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A phase-field approach describing the dynamics of a strained sohd in contact with its melt is devel- 
oped. Using a formulation that is independent of the state of reference chosen for the displacement 
field, we write down the elastic energy in an unambiguous fashion, thus obtaining an entire class 
of models. According to the choice of reference state, the particular model emerging from this 
class will become equivalent to one of the two independently constructed models, on which brief 
accounts have been given recently [J. Miiller and M. Grant, Phys. Rev. Lett. 82, 1736 (1999); K. 
Kassner and C. Misbah, Europhys. Lett. 46, 217 (1999)]. We show that our phase-field approach 
recovers the sharp-interface limit corresponding to the continuum model equations describing the 
Asaro-Tiller-Grinfeld instability. Moreover, we use our model to derive hitherto unknown sharp- 
interface equations for a situation including a field of body forces. The numerical utility of the 
phase-field approach is demonstrated by reproducing some known results and by comparison with 
a sharp-interface simulation. We then proceed to investigate the dynamics of extended systems 
within the phase-field model which contains an inherent lower length cutoff, thus avoiding cusp 
singularities. It is found that a periodic array of grooves generically evolves into a superstructure 
which arises from a series of imperfect period doublings. For wavenumbers close to the fastest- 
growing mode of the linear instability, the first period doubling can be obtained analytically. Both 
the dynamics of an initially periodic array and a random initial structure can be described as a 
coarsening process with winning grooves temporarily accelerating whereas losing ones decelerate 
and even reverse their direction of motion. In the absence of gravity, the end state of a laterally 
finite system is a single groove growing at constant velocity, as long as no secondary instabilities 
arise (that we have not been able to see with our code). With gravity, several grooves are possible, 
all of which are bound to stop eventually. A laterally infinite system approaches a scaling state in 
the absence of gravity and probably with gravity, too. 

PACS numbers: 81.10.Aj, 05.70.Ln, 81.40.Jj, 81.30.Fb 
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I. INTRODUCTION 



Already when introducing the notion of a surface quantity Gibbs imphcitly entertained the idea 
of a phase field 0: any density of an extensive quantity (e.g., the mass density) between two 
coexisting phases changes gradually (but swiftly) from its value in one phase to its value in the 
other. The existence of a transition zone, though microscopically of atomic extent (far enough from 
a second-order phase transition), underlies the very Gibbs definition of surface quantities. In phase 
transition phenomena, either of first or second order, this notion has been adopted in Landau's 
spirit. Because energy is an extensive quantity, too, there is an extra energetic cost associated 
with the transition region, characterized in the appropriate thermodynamical potential density by 
a term of the form e*(V(/))^, e* being the stiffness of the transition region. 

The notion of a phase field has appeared abundantly in the literature in the context of phase 
transition phenomena The transition width diverges for a second-order phase transition 

at the critical point, and thus it is essential to introduce the transition region. For a first-order 
transition, such as a liquid-solid interface, confering an importance to the interface thickness may 
seem quite anecdotic if one is interested in properties which occur on a scale larger than the atomic 
one; typical examples are dendritic patterns occuring at the scale of a /im. Nevertheless, it is here, 
where phase-field modeling has become most useful in numerical treatments. 

Before phase-field models became popular, it seemed quite natural to treat the surface as a 
geometric location on which boundary conditions are imposed (e.g., for a moving front the normal 
velocity is proportional to the jump in the gradient of the temperature or concentration field). 
This is the so-called sharp interface approach, adopted both in analytical and numerical studies in 
a variety of contexts of front problems. 

There has been an upsurge of interest in the phase-field approach to free-boundary problems 
more recently, though the method was actually introduced pretty early |^,^, as a computational 
tool to model solidification. Various studies |qj7|,p|,p|,p^ 11 1 have demonstrated the virtues of this 
method in moving-boundary problems. 

Regarding the way how to use phase-field models, there are two distinctly different philosophies. 
These may be best discussed considering dendritic growth, where a set of well-established contin- 
uum equations exists describing phenomena in terms of a sharp interface. On the basis of this 
knowledge, a phase-field model can be justified by simply showing that it is asymptotic to the 
correct sharp-interface description, i.e., that the latter arises as the sharp-interface limit of the 
phase-field model when the interface width is taken to zero. This is definitely a sufficient condi- 
tion for the phase-field model to yield a correct description of the continuum limit, providing the 
interface thickness is taken small enough. Small enough sometimes may mean impractically small. 
The second approach to phase-field modeling is to guess or derive an appropriate form for the free 
energy of the two-phase system, including the energy cost of the transition region and to regard 
this as a physical model in its own right. In this case, one might actually forgo considering the 
limit of small interface thickness and such a model would even make sense, if the strict limit of 
vanishing interface width did not correspond to sensible physics. 

Of course, a problem arises, if a phase-field model obtained in the second way gives predictions 
that are different from that of the sharp- interface equations. A follower of the first philosophy would 
then discard the phase-field model, whereas one of the second might contemplate the possibility 
that his model contain more physics than the sharp-interface model. In the case of dendritic 
growth the situation is pretty clear: the sharp interface model gives the right answers. However, 
this statement cannot be generalized easily, since not all sharp-interface models are as well-founded 
as that for dendritic growth and because the extreme smallness of the interface width cannot be 
always guaranteed (it might for example become doubtful for a phase transition that is only weakly 
first order). 

A related issue is the question of thermodynamic consistency, i.e., the derivation of the model in 
the spirit of Gibbs from a free-energy or entropy functional. It is clear that with a known sharp- 
interface limit in mind, there is no need at all to obtain a phase-field model this way (which would 
mean to make it "thermodynamically consistent" ) as long as one ensures its asymptotic approach 
to this limit. In fact, it has turned out that in some cases where both a thermodynamically 
consistent formulation of a phase-field model and a nonvariational formulation exist, the latter was 
numerically more efficient [n2| and hence preferable on practical grounds. 
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On the other hand, thermodynamic consistency has its virtues. This can be seen particularly 
well in the case considered here, the influence of elasticity on the stability of a solid interface. It is 
quite straightforward to write down the contribution of the elastic energy to the total free energy. 
Hence, if we have a good idea about the physical origin of the free energy to be considered, the 
corresponding phase- field model is easily obtained, and it is bound to be right. As a result, one 
may derive sharp-interface equations in cases where they are not known. 

For the Grinfeld instability to be considered here, the sharp-interface equations are well known. 
Nevertheless, it is of course tremendously satisfying, if they simply pop out of the phase-field 
equations as the sharp- interface limit. Not only does this provide a natural countercheck of our 
ansatz for the free energy, but it also gives us a new angle of view at the instability, leading to the 
prediction of circumstances in which the Grinfeld instability should not occur under anisotropic 
stress, but might appear with isotropic stress. We shall consider this point in section II. 

Let us return to the advantages of the phase-field method. The first virtue of phase fields is 
pretty obvious: instead of tracking permanently the a priori unknown interface position in the 
sharp-interface limit, and imposing nontrivial boundary conditions for the discontinuity of the 
fields, the interface in the phase-field approach is nothing but the location of a rapid variation 
of the field 0, while the two phases are treated as the same entity. Thus there is no boundary 
condition to be imposed in the transition region, a fact which greatly facilitates both the numerical 
implementation and analysis. This is done at some price: one must, in principle, mix disparate 
length (and thus time) scales: the pattern length and the interfacial width, whose ratio may range 
over many orders of magnitude. This may render the numerical procedure excessively expensive, 
a fact which would quickly take us back to the sharp-interface problem, where the small length 
scales out of the problem. 

As discussed recently [ p3[ , the sharp- interface limit that one would like to represent when writing 
the phase-field equations, only makes physical sense on "outer" length scales much larger than the 
physical extent Ic of the transition region, and thus does not depend structurally on the details of 
the interface shape on the inner Ic scale. The mathematical question, formulated in the framework 
of phase-field models, of formally recovering the sharp interface description via an asymptotic 
(multi-scale) expansion in the limit Ic might, from this point of view, seem irrelevant. 

It should however be kept in mind that the actual matching conditions arc imposed for the 
limit, where an "inner" variable (defined in the transition region) goes to infinity. This entails a 
certain amount of liberty in the choice of functions defining the free-energy density, because the 
precise behavior of these functions on the inner scale does not matter. Hence, the validity of a 
phase-field model can indeed be judged by simply showing that it asymptotically reproduces the 
correct sharp-interface description. Whether the additional information encoded in the structure 
of the phase-field model on the inner scale is physically relevant, is a question to be decided on a 
case by case basis (as implied by our discussion above). 

An example where this is relevant, is provided by the Young condition for the contact angle 
of a droplet on a substrate. This is a condition on "outer" scales, while the inner scale is rather 
governed by van der Waals interactions of a thin liquid film with the substrate, leading to some 
nontrivial corrections of the wetting profile at small atomic length scales. In other words, what 
matters in the phase-field description is not that the width of the interface be of atomic extent, 
but rather that it be small in comparison with the scale of the pattern of interest. 

This physical argument has been cast in a mathematical form in [ p^ , where the thin-interface 
limit was considered, arising from an alternative asymptotic procedure. This has to be contrasted 
with the sharp-interface limit with the small parameter being the ratio of the interface thickness 
to the capillary length. Making the width of the interface small only in comparison with the scale 
of the pattern leads to a rather important enhancement of the computing speed, thus rendering 
the phase-field approach attractive with regard to numerical efficiency as well. Unfortunately, 
in most systems the thin-interface limit is not as easily accessible as for dendritic growth in the 
thermal model [ p^ . Hence, it is difficult in general to make our argument mathematically rigorous. 
However, we will give a consideration of length scales in section II that clarifies what are the ratios 
that have to be kept small for our model to be a good description. Moreover, we have checked 
that the dependence of the results on the interface width becomes weak for the small values that 
we use for the latter. 

Finally the phase-field approach has the additional virtue of regularizing instabilities, such as 
the development of cuspy structures, often setting severe limitations in numerical studies. In the 
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elastic system, we are led to the question whether the sharp-interface models still make sense in 
the limit where they predict finite-time singularities. No such singularities arise in the phase-field 
model, thus possibly extending the range of validity of the latter beyond that of the former. This 
will be discussed in more detail in section IV. 

In the context of growth phenomena phase-field approaches have been introduced in problems 
involving temperature or concentration fields. There are myriads of situations, however, where the 
corresponding transition is monitored, or at least affected, by strain. A typical situation is a solid 
under uniaxial stress. This leads to the Asaro-Tiller-Grinfeld (ATG) instability (see Ref. 

| |l6[ | for a recent review. A surface corrugation allows to lower the stored elastic energy. Other 
examples of particular interest include solid-solid transformations, phases in nonequilibrium gels, 
molecular beam epitaxy, solidification of lava, etc. It is thus highly desirable to develop a phase- 
field approach including the stress as an active variable. Very recently, two groups, Miiller- Grant 
(MG) and Kassner-Misbah (KM) (l^Jl^ , have independently developed such an approach and have 
given a brief account on it. 

This paper will present extensive discussions of this question and give new results. We shall also 
provide a comparison between the MG and KM models, and point out similarities and dissimi- 
larities. As already known from sharp-interface simulations a solid under stress presents 
a somehow stringent behavior in that no stable steady-state solution seems to exist. This is also 
the case from analytical studies in the long- wavelength limit Nozieres had shown that 
the bifurcation from the planar front to the deformed one is subcritical (the analogue of a first- 
order transition). The study of Spencer and Meiron [^3|] focused on structures with a given basic 
wave number in the absence of gravity (where the instability does not have a threshold) and on 
systems, in which the transport mechanism necessary for the instability to manifest itself is sur- 
face diffusion. They find that in the unstable range of wavenumbers (i.e., for wavenumbers below 
the marginal value, above which surface tension stabilizes the planar interface) there exist finite- 
amplitude steady-state solutions, if the wavenumber is close enough to marginal. This branch of 
steady-state solutions terminates by structures developing cusp singularities, despite the stabilizing 
influence of surface tension. It cannot be overemphasized that this result is not an artifact of their 
numerics. Indeed, they investigate carefully the effects of numerical fine-graining using a code with 
spectral accuracy, and their discretization sequence seems to get fine enough to render their ex- 
trapolation valid. The evidence for the appearance of true cusps in the sharp-interface continuum 
model becomes compelling, if one takes into account the work of Chiu and Gao |24| who found 
analytical solutions developing cusp singularities in finite time. The conclusion of Spencer and 
Meiron is that for generic initial conditions, including sufhciently small wavenumbers, finite-time 
singularities will always occur. Moreover, they state that this is also true in the presence of gravity 
(beyond the threshold) and under fairly general conditions. 

For a physical system, finite-time singularities will be prevented by intervening effects that are 
not considered in the model. This would mean that nonlinear elasticity and plasticity have to be 
taken into account. For example, the formation of dislocations (plasticity) could blunt cusps again. 

Two questions then arise. What kind of structures can be expected before cusp formation and 
what kind of structures will prevail eventually? Our previous study has shown that the initial 
cellular pattern may develop into a super-structure, where a groove or several of them accelerate 
in a spectacular fashion, thus relieving the stresses in their surroundings significantly and causing 
nearby grooves to recede. Further evolution of the structure was difficult to handle numerically, 
due to the development of cusps which appear in the numerics as they should according to the 
analytic results In the phase-field approach presented here, no cusps can arise. The question 
is of course legitimate, whether the model, which allows to track the dynamics of the structure 
much beyond the times where earlier studies had to fail numerically, still gives a faithful description 
of the physics. Here we take the point of view, that the details of the description in the locations, 
where stresses become large, may not be correctly captured by the model, since we neither include 
nonlinear elastic effects nor other effects such as capillary overpressure explicitly ||2^ . If they have 
the right sign, these might prevent cusp formation (23|. However, the result of any such effect must 
be to blunt cusps, which the phase-field model does. As we shall see, it does so in a non-obtrusive 
way by introducing a cutoff for interface curvature. Moreover, away from the cusps, stresses are 
low enough for linear elasticity to apply. Hence, we believe that the development of the overall 
morphology is still correctly described by the phase-field model. 

A more detailed justification would point out that nonlinear elasticity will first make itself felt 
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by stresses increasing more slowly as a function of strains than in the linear case and that next 
plasticity will act to introduce an upper cutoff for stresses, where the material will yield. Now 
the effect of any resulting modification in the stress-strain relationship on the remaining body 
can be reproduced by cutting out the piece of material where linear elasticity ceases to hold and 
by requiring boundary conditions at the edge of the cut-out piece that correspond to the correct 
stresses. In the case of material yielding near a singularity of the stress field (obtained assuming 
linear elasticity), these boundary conditions will essentially be that the stresses are close to the 
yield stress at the boundary. This is mimicked by the phase-field model in which the maximum 
supportable stress is, for a given geometry, determined by the interface width. 

Therefore, we think this is one case, where the phase-field model can do more in the description 
of the physical system than its sharp- interface limit. 

It will emerge that usually one leading groove continues to deepen while neighboring ones recede 
after the winner has started to relieve the stress that kept them growing. Finally the surface shows 
a single deep groove evolving in time and becoming a location of a strong stress accumulation, 
possibly until the fracture threshold is reached. Presumably, before that stage is reached the 
validity of the model will break down. We shall make some speculation on future directions to 
elucidate the physical behavior in real systems. For generic initial conditions, we may consider the 
dynamics a continual coarsening process which initially develops as described in and later is 
dominated by groove growth and shrinkage. 

The paper is organized as follows. In section II we give the continuum equations ordinarily 
used in the description of the Grinfeld instability. This is mainly done in order to introduce 
appropriate length and time scales in nondimensionalizing the equations; then we present our 
phase-field approach and discuss how the interface width has to be chosen in comparison with the 
other length scales. We demonstrate how the phase-field model can be employed to derive new 
sharp-interface equations in the presence of body forces breaking rotational invariance. 

Section III presents validation results, a comparison of the MG and KM models and describes the 
main findings of our simulations. Section IV sums up the results and discusses perspectives. The 
mathematically rigorous asymptotic expansion used to derive the sharp-interface limit has been 
relegated to an appendix, as the calculation is somewhat lengthy and would interrupt the fiow of 
the text. Since we use the MG model in a slightly different form from that presented originally 
jirt , we give the connection between the two formulations in a second appendix. Finally, appendix 
C contains the analytic derivation via conformal mapping of the stresses for a particular interface 
shape to be compared with the numerics. 



II. THE GRINFELD INSTABILITY 



A. Sharp-interface equations 

A description of the basic ingredients of the Grinfeld instability has been given elsewhere . 
Therefore, we may restrict ourselves to explaining the physical mechanism and giving the equations. 

We wish to describe the behavior of a solid submitted to uniaxial stress, at the surface of which 
material transport is possible. Consider the example of a solid in contact with its melt. An 
accidental corrugation of the surface will act to reduce the stress at its tip and increase it in the 
valleys next to it. That is, the solid can decrease its elastic energy density by growing tips (where 
the stress is lower) and by increasing the depths of valleys (where it gets rid of material having 
a higher density of elastic energy due to larger stresses). This tendency is most easily cast into 
equations by writing down the chemical potential difference A/i = /is — ui at the interface (the 
subscripts refer to the solid and liquid or nonsolid phases, respecively) 

A/i = — — [att - cr„„) H 7K -f gCix) (1) 

^i^Ps Ps Ps 

Herein, the first term is of elastic origin; ctu and cr„„ are the normal stresses tangential and 
perpendicular to the interface, E is Young's modulus, v the Poisson number, and ps the density 
of the solid. The second term describes the stabilizing infiuence of the surface stiffness 7, taken 
isotropic here (so it becomes identical to the surface energy), k is the curvature of the interface; for 
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simplicity, we consider the two-dimensional case only. Finally, the third term is the contribution of 
gravity (g), where Ap — Ps — pe is the density contrast between the solid and the liquid (or vacuum) 
and (^(cc) is the interface position, given by its z coordinate (the z axis is oriented antiparallel to 
the gravitational force). Equation (|^) holds for plane strain. For plane stress, the prefactor 1 — 
has to be dropped. 

The dynamics is then described by giving the normal velocity in terms of the chemical potential 
difference. For a solid in contact with its melt this would simply be 

Vn = -^Afl, (2) 

where k is an inverse mobility with the dimension of a velocity. In the case of a solid in contact with 
vacuum and surface diffusion as the prevailing transport mechanism we would have Vn — D'S/^Ap 
instead. 

Of course, in order to compute u„, we must first obtain the stresses entering (Q). This involves 
solving an elastic problem (djaij = 0) with a prescribed external stress and boundary conditions 
on the interface, assuming an appropriate constitutive law. Ordinarily, Hooke's law for isotropic 
elastic bodies is assumed [therefore we have only two elastic constants in (|l|)]. Neglecting the 
capillary overpressure, which usually is a good approximation, we have as boundary conditions at 
the interface (T„„ — ~p, where p is the pressure in the second phase, and ant = 0, i.e., the shear 
stress vanishes. 

A linear stability analysis of a planar interface under the dynamics given by equations (|l|) and 
(H) yields the following dispersion relation {oj is the growth rate, q the wave number) 

1 f 2^2(1 



c. = ^|^V-^.-7.^-Ap.9r (3) 



(To is the uniaxial external stress. Eq. (^ provides us with a critical wave number Qc = \J Apg/7 

(an inverse capillary length) and a critical stress ctqc = \^fqcE j (1 — j/^)] ^ , below which the planar 
interface is stable. The wave number of the fastest-growing mode can be inverted to give a length 

. - (4) 



Apart from a prefactor of 2/7r, this length is identical to the so-called Griffith length. 

Note that a planar interface will not remain at its original equilibrium position, even when only 
a subthreshold stress is applied, i.e., when it is "stable". Our dynamical equations predict that it 
has a nonzero velocity as long as A/i is different from zero. Hence it will recede to smaller values 
of If the density of the solid is bigger than that of the liquid, the chemical potential on the 
solid side of the receding interface decreases faster than that on the liquid side and there is a new 
equilibrium position, which can be computed directly from (|l|) and which evaluates to 

Equations (^ and (|^) provide us with two independent length scales of the problem, the first 
of which is due to a competition between stress and surface energy, while the second arises from 
the competition between stress and gravity. For the purpose of nondimensionalizing equations, t\ 
is more appropriate, as this length does not diverge in the limit of vanishing density difference 
(or gravity). To obtain a natural time scale t, we can replace q in either of the two wave- number 
dependent terms of (0) by l/^i. This leads to 



(6) 



The nondimensional version of the dispersion relation then reads (a) = rw, q = £iq): 

^^2,-,--!-. (7) 
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which shows clearly that the problem without gravity (when £2 becomes infinite) can be made 
parameter free, i.e., elastic and other parameters only set the time and length scales; apart from 
that we should expect the same dynamics for all systems. With gravity, the dynamics is essentially 
determined by the ratio of the two length scales introduced. 



B. Elastic energy and state of reference 

Let us now proceed to investigate the contributions to the free energy of the same system. The 
phase-field model will then consist in writing the free-energy density that takes into consideration 
the global elastic energy in both phases. 

As is usual with elastic problems, it is important to specify the state of reference defining 
the positions of material particles with respect to which displacements are measured. This is 
crucial whenever the reference state is not that of an undeformed body but one that is subject to 
prestraining (which will turn out useful later). In order to make this point clear, and in the hope 
of helping subsequent discussions, we would like first to dwell on this issue. 

If the only energy present is elastic energy and Hooke's law holds true, the free energy per unit 
volume can be written as 

/ = /iUy Uy + , (8) 

where summation over double subscripts is implied. A and fj, are the Lame constants, fi being 
better known as the shear modulus. For plane strain, these elastic constants are related to Young's 
modulus and the Poisson ratio via /i = E/[2{1 + 1/)] and A = Ev/[{1 + v){l — ^v)]. The stress 
tensor aij is then 



d f S ' 

^— = 2nuij + XukkSij = 2/x(u,y - Ukk-j) + KukkSij . (9) 



K — 2/i 4- A/d is the bulk modulus, and the last relation has the advantage of making explicit 
the parts of the stress tensor causing pure shape and pure volume changes, respectively, {d is the 
spatial dimension.) We will nevertheless mainly use the relations containing ^ and A, which are 
more compact. The implied reference state here is Uij = 0, for which aij = 0. 

However, if we choose a reference state given by a different strain tensor setting u.y — 

Uij — uf^^ , then we should not simply replace Uij by iiij in (|9|) , as the stress tensor and the elastic 
energy are, in principle, measurable quantities and should thus be unaffected by a change of strain 
reference state. Hence, we have to write 

<7y = 2^l {u,j + uf^^ + A {ukk + wife ) 5^j , (10) 

and change definition (^) accordingly, i.e. replace Uij by Uij +uf^ . In this situation, the zero-strain 
state would not be stress-free. An alternative way to specify a reference state would then consist 
in giving the stress of the zero-strain state. 

In general, the thermodynamic system under consideration will not only contain elastic contri- 
butions. Then the equilibrium state, corresponding to a minimum of the free energy, may not be 
a state of vanishing strain. A trivial example is a solid in equilibrium with its melt, where the 
equilibrium state in the solid corresponds to the strain produced by the equilibrium pressure p of 
the liquid (the equilibrium stress tensor of the solid is —pSij). The form of the free-energy density 

accounting for such a situation is not (||) [which does not exhibit a minimum at u^-"'^^ 0] but 

/ = ^^ - .g'^') (u., - u^^) + \ [uu - v^)" . (11) 

This is manifestly minimum at u^i^^ and the nonzero value of the latter quantity takes into account 
nonelastic contributions to the free energy. If we now define the stress from the first relation in 
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Eq. (^), i.e., aij — df /duij, it will be nonzero only, if there are forces driving the system away 
from equilibrium. If we rather define it via the second equality of Eq. (^), meaning that we set 
aij — 2iiUij + Xui^kSij (which is now different from df /duij), it will describe, in addition to these 
forces the prestress necessary to keep the system in equilibrium. An invariant relation between 
stresses and strains follows from requiring 



/=/._, K-^S'Mrf«.., (12) 



which leads to 



2m (u« - u,'r') + - «£") 4i ■ (13) 



Once both cr^^'^' and u''^'^^ are specified, this equation gives us an unambiguous relationship between 
stresses and strains. Depending on which variables we choose to define the reference state, we obtain 
the conjugate variables of the same state from (|l^). If we choose, for instance, a strain- free state 
as reference, this equation will provide us with the corresponding stress of reference, if we choose 
a stress-free state of reference, it will yield the strain of reference. 

As an example we can look at a case where the equilibrium stress is cr|°'''' = —poSij, and ask 
what should the strain be. Hooke's law - written in a such a way that the absence of strain implies 
the absence of stress as well [see @] - then gives us an equilibrium strain 

Since the free energy must not depend on the choice of reference state, it is clear that it does 
not matter whether we use Uij or Uij in providing that we use the correct values of uf^^ 

and uip'' , respectively. Suppose that we choose another reference state characterized by the strain 
tensor Uij in such a way that when the strain is zero, the stress is equal to —posSij. A vanishing 
strain then corresponds to a pre-stressed situation. If the equilibrium stress is again —poSij as 

above, we must have a new equilibrium strain u[j'^^ obeying, according to our invariant relation 
(p^), —posSij + PoS'ij — —2^uf^^ — Xu'^f^^Sij. After a simple manipulation we obtain 



-(cq) _ POs_PO^ 

Of course, Eq. ( |l^ ) is a special case of ( [l5| ) for pos — 0. The free energy, expressed by Uy , then 
reads 



= ^iUijUij + ^UiSn + (Po - PQs)uii + 2(2^ + \d) ~ ^^'^"^ ^^^^ 

It should be realized that Eqs. ( |ll| ) and ( p^ describe exactly the same situation, if corresponding 
values for the strain fields without and with a tilde are inserted. At this point we have said 
nothing about applied external stresses or so. However, if we choose, say, vanishing displacement 
as boundary condition at a planar interface, directed along the x direction, then this will correspond 
to two different physical situations for the two different choices of the state of reference. Let us for 



simplicity assume pq — 0. According to (|14D, we then have u^"^'' = 0, and (|l l| ) implies (^. Setting 



Uxx = in Eq. (bf) we obtain, because of the boundary condition Uzz = that also u^z — 0, and 
there is no stress at all. On the other hand, if we set u^x = 0, we have to use (p^ ) and (p^ ) to obtain 
the elastic state of the solid. The boundary condition for Ozz implies Uzz — Pos/i'2^J■ + A), which 
in turn leads to axx = — 2/ipos/(2/i -f A); hence vanishing displacement along our planar interface 
means a solid that is homogeneously strained in the x direction with a prestress (Tq = f^xx- 

As we shall see later, the latter choice of the state of reference has been made in the phase-field 
model discussed in [0, the former (setting pos — 0) in [Q. These are the most natural choices, 
although an infinity of (less natural) alternatives is available. 
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C. The phase-field model 



The total (solid+liquid) free energy of the system can be written as 



F[^, {Uy}] = / dV 



(17) 



where e is a length parameter controlling the order of magnitude of the transition region described 
by the phase field. T = 37/e is the energy density corresponding to the surface energy 7 being 
distributed over a layer of width « e. (The factor 3 is just a convenient choice, simplifying later 
derivations.) 

If we start from the invariant form (|l^), we can set up a whole class of phase- field models 
at once and specify the reference state later. In order to be able to write a single elastic energy 
expression for the two-phase system, we formally treat the liquid as a shear-free solid (not including 
hydrodynamics). We will discuss some consequences of this approach later. 

A straightforward ansatz for the elastic energy density is then 

/el(0, {Uy }) = h{(j))fso\{{u^j}) + [1 - /l(0)]/liq({"ij }) 

+ , (18) 

where fso\{{uij}) and /iiq({'Uij}) are the densities of elastic energy in the solid and in the liquid, 
respectively, and where h{(j)) may be interpreted as a "solid fraction", which must be equal to one 
in the solid and equal to zero in the liquid. We choose h{<j)) = (f)^ {3 — 2(j)) for reasons of convenience: 
with this choice h'{(f)) — for </> = and for 0=1, i.e., in the bulk phases. This leads to the 
advantage (see appendix A) that the zeroth-order solution of the asymptotic expansion in powers 
of e is valid to all orders in the outer region considered. 

Since different reference states may be chosen in the solid and in the liquid, the equilibrium 
strains carry a subscript s or £, respectively. This would not be necessary here, because the 
prefactor [h((f>) or 1 — /i(0)] decides whether the equilibrium expression for the strain in the liquid 
or in the solid has to be taken. However, as soon as we take derivatives with respect to 0, this 
criterion of distinction becomes ambiguous, so we prefer to make the difference explicit from the 
outset. A is the bulk modulus of the liquid. 

To account for the possibility of a phase transition, we introduce a double well potential 

/dw(0) = 2rg(0) , (19) 

where g{(j)) — (?!)^(1 — (p)'^. The minimum at </> = 1 corresponds to the solid phase, the one at (/> = 
to the liquid phase. Note that while this potential looks similar to the fourth-order polynomial 
used in the Landau theory of second-order phase transitions, it is employed in quite a different 
manner here. The two minima correspond to the two phases and the symmetry of the potential is 
of secondary importance; in Landau's approach, symmetry considerations are at the heart of the 
theory, the symmetric minima describe the same phase, and the second phase corresponds to the 
unstable maximum in between. Since in our case both phases sit at a minimum, the transition 
described by the double well potential is of first order. We do not need a sixth-order polynomial 
as would be necessary in Landau's theory for first-order phase transitions. 

Gravity will be included in essentially the same way as in the sharp- interface equations discussed 
above; i.e., its effect as a body force in the mechanical equilibrium condition is neglected but its 
influence on the chemical potential is taken into account. This is a good approximation usually (one 
can estimate the cross-effect of gravity on the elastic energy to be on the order of psgH/ao ^ 1, 
for typical heights H of the sample) . Then the contribution of gravity to the free-energy density 
becomes 

/grav(0, ^) = (z - ZQ){psh{(j)) + pe [l~h{4>)]}g 

= {z- zo)h{(j))Apg + {z- zo)ptg , (20) 
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where we have taken the zero point of this potential energy at z = zq. Note that in taking fixed 
values for Ap and for ps, we also neglect the second-order effect caused by density changes due to 
strain. 

Finally, we wish to be able to control the equilibrium position of the interface "by hand" via 
addition of a constant to the free-energy density of one phase; this phenomenogical contribution 
to the total free-energy density may be conveniently written as 



= -hi<t>) 



1 



2E 



2p + X 2 
8p{p + xf'° 



(21) 



and it is normalized such that setting <tqq = ctq , we can keep the equilibrium position of the planar 
interface at the fixed value zq, independent of ao- This is useful, for example, if one wishes to assess 
the relative position of the maxima or minima of an evolving structure with respect to a planar 
interface at the same external stress. Because of the recession of a planar interface according to 
(j^), such a comparison would otherwise be difficult. 

Collecting all contributions, we obtain for the total free-energy density 

f{(j> , {Uij}, Z) = /d„('?!)) + /ol((/', {Uy }) -I- /grav(</', z) -f /c(0) 



+ [1 - hm \ {y^u 



,(cq) 



Apg{z - zo) 



2^ + A 
8Ai(Ai + A)' 



+ (z - ^o) 



(22) 



Note that here the terms in braces, in particular the elastic term, have acquired a prefactor e. This 
e dependence is spurious, as we have taken the prefactor F oc 1/e in front of everything, and the 
factor e just serves to cancel this. In fact, the only contribution to the free energy that can depend 
on e explicitly is the double well potential, which must ensure that in the limit e ^ the only 
possible states are the bulk phases and must therefore become infinite for all values of different 
from 1 or 0. All the other energies can depend on e only implicitly via /i((^), the local solid fraction 
of the two-phase system. 

We then require </> to satisfy a relaxation equation for a non-conserved order parameter. This 
equation takes the form 



dt 



-R- 



.6F 



(23) 



and the prefactor R should contain the mobility 1/fc defined in (||). The dimension of R must be 
(energy density x time)~^, which leads us to choosing R = l/(3fcpse). This essentially amounts to 
setting the time scale for the evolution of (j) (which must be related to the width of the transition 
region, because it is only in this region where has an appreciable dynamics). 
We arrive at 



dt 



7 

kps 



-l(2.g'(0) + ^;i'(</.){^( 

( 
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A 

2 



S)(' 



,(oq) 



[Uii ~ I 

Jp+^ ,\ 
8M/i + A)''™ij 



(cq) 

ii.e 



+ iz- zo)Apg - 



(cq) 



(24) 



Herein, g'{(f>) and h'{(f)) are the derivatives of /i('/>) and g('/>) with respect to their argument . As we 
have mentioned before, h'{(p) vanishes in the solid as well as in the liquid phases [see eq. (A24)]. 

In writing down an equation for the evolution of the elastic variables, we have to be careful 
about the fact that the strains Uij, i,j — l,...,d are not independent quantities. Therefore, the 
variational derivatives SF/Suij are not independent. Instead of introducing Lagrangian multipliers, 
we can however exploit the fact that the components Mi, i = l,...,ci of the displacement, related 
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with the strains via Eq. (|33|), are independent variables. Assuming that the time scales of our 
problem are large in comparison with sound propagation times, wc conclude that the variational 
derivatives SF/Sui are equal to zero. This is is an adiabaticity assumption. Hence, we obtain 

Q_ SF d 6F 



5ui dxj Suij 

= {hW {<^^, - 4"^) - [1 - hm [p - p^^'^^) % } . (25) 

This is nothing but a generalized elasticity problem, with the generalized stress tensor given by 
the quantity in braces. 

Before moving to a demonstration of the sharp-interface limit, let us discuss scales. Since the 
elastic problem ( p5| ) is formally linear in the strains, rendering it nondimensional is straightforward 
and unenlightening. On the other hand, trying to cast (p^) int o non dimensional form, we realize 



that besides the length and time scales discussed in subsection |II A| , we need a third length scale 
£^ = "f/K, apart from the width of the transition region e. So the phase- field model contains four 
length scales altogether. 

Normalizing elastic moduli and stresses by the bulk modulus, i.e., setting M =^ fi/ K, A = \/ K, 
A — X/K, Soo = o'oo/K, we obtain 



dt 



- § (2.9'(0) + 3^/^'(0){m - 4-)) - 4-)) + I {uu - )' 
A/ (cq)\2 24.. ^M + A_ 2 \A 



(26) 



where t — t/r is the nondimensional time, and z — z/£i, V — iiV are nondimensionalized spatial 
operators. Physically, represents an atomic scale. For many materials, "f/K is on the order of 
the lattice constant. 

For the phase-field model to work properly, we must impose some conditions on the length scale e. 
We definitely need e/£i <C 1 to have a decently sharp interface. Moreover, the /i'(0) term must not 
become too large in comparison with the g'{(t)) one, otherwise the minima of the double well move 
away from the positions (j) — and (f> — 1. This appears to suggest that we also need e/i^ ^ 1. We 
compute some typical values. Using the material parameters of solid He |2^], the system for which 
the Grinfeld instability has been unambiguously demonstrated by Torii and Balibar ]2^ ] , we obtain 
the estimates £i « 0.1 cm, £2 « 0.1 cm, £3 « 10~^ cm, and t w 1 s. If we had to require e <C 4, we 
would have a problem with very disparate length scales, as our numerical grid would have to be 
smaller than £3, whereas the length scales governing pattern formation are £1 and £2- Fortunately, 
the quantity e/3^3 appearing in (|2^ ) is multiplied by squared strains, and the are on the order 
of 10^^. Moreover, we have 2£3/£2 ~ 2 x lO^'* and the last term in braces can be estimated by 
|;Sqq w 2 X lO^'*. Therefore, the actual condition for our model to be useful is 10^* x e <C 8^3, 
which is much easier to achieve. In our simulations, we typically had 10~^ x e/3^3 « 0.1. 



Equations (|2J)-(25) constitute the basic phase-field equations for the phase transformation under 
stress. 

To specify our model completely, we have to indicate the equilibrium stresses and strains. Let 
us assume the following forms for the stress-strain relationships in the two phases, 

CTij = -posSij + ^fiu^j + XukkSij , (27) 
P^Poe- >^Ukk , (28) 

and require the equilibrium pressure to be po- For a planar interface, this fixes the normal stress 
in z direction to be cri^'^^ = —po- If we assume the equilibrium stress tensor to be isotropic (a very 
natural assumption in most cases), we have cr^^'^^ — —poSij, and in the solid is given by (p^. 

In the liquid, we have u'^^^'f = {po£ — po)/^- Note that only the displacement divergence Vu = ua 
appears in the elastic energy of the liquid. This gives us a degree of freedom (neither Uxx nor u^z 
are fixed separately in the liquid, only their sum is) that will turn out important later. (Without 
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this degree of freedom, it would not be feasible to treat the liquid as a shear-free solid, as will be 
discussed in appendix A.) 

Inserting these equilibrium values into (|2j) and (p5|), we obtain as basic equation of motion for 
the phase field (introducing the abbreviation k = kps) 
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(29) 



where we have defined further abbreviations 

Ap = poe - pos (30) 

and 

AW ^^(PO-POsf 1 (PO - POi? 2fl + X 2 (ou 



,2 



The elastic problem can be cast into the suggestive form 

dx 



^^4r{ K4>)^r3 - [1 - h{4>)]p5^, ) , (32) 



from which it is even more transparent that the expression in braces is nothing but a generalized 
stress tensor. Note that the phase-field model always guarantees exact mechanical equilibrium 
with respect to this stress tensor, but that the validity of a linear relationship between strains and 
generalized stresses is only warranted outside the interface region, where the values of (j) cease to 
depend on the Uij. (This means that in the vicinity of sharp groove tips we will automatically 
have deviations from Hooke's law, albeit they are not modeled to satisfy a particular nonlinear 
constitutive relation.) 

These equations arc to be solved subject to the conditions that the phase field approaches its 
limiting values in the bulk phases. To make them closed equations, we have to replace (7^ and Uij 
by the field variables Ui using the definition of the strain tensor, 

1 / du, duA 

and Hooke's law. 

It remains now to be shown that this model reproduces the sharp-interface limit when the width 
of the interface is small. This calculation is given in the appendix. Its central result is formula 



(A50), which we rewrite here in nondimensional form (for ctoq — 0) 

1 {att - Onnf , - , ^1 



2 al '-24(^-^°)i- 

The ansatz proposed in |Q is slightly different. One difference mainly concerns the interpreta- 
tion or philosophy of the approach. In the MG model, the phase field is considered the variable 
determining the shear modulus. The shear modulus is the macroscopic quantity deciding whether 
a piece of condensed matter is solid or fluid. Hence, the phase-field order parameter differentiates 
between liquid and fluid and has a transparent meaning in the context of liquid-solid transitions. 
Of course, the model can be extended easily to the case of two solids with nonvanishing shear 
moduli on both sides of the interface. In the KM approach, the traditional and more conventional 
view is taken that the phase field decides between two phases characterized by their respective 
free energy densities. That one of these phases is a liquid is of secondary importance, as it were. 
Again, in principle, it might be another solid. Of course, i/ the second phase is chosen a liquid, 
then its shear modulus must vanish. And indeed, this is guaranteed in the current form of both 
models by construction. For ease of further comparison, we give the phase-field equations of |0 
in appendix B and show how they are mapped onto the form ([29|j3^). 
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In concluding this section, we would like to comment briefly on the consequences of an anisotropic 
equilibrium strain. Suppose we submit a body consisting of piezoelectric material to a homogeneous 
electric field. (Alternatively, we could consider some magnetrostrictive material under the influence 
of a magnetic field.) This body will contract or expand until it reaches a new equilibrium state 
compatible with the body forces exerted by the field. The new state will have anisotropic strain 
and, assuming isotropic elastic properties, an anisotropic stress tensor as well What will the 
surface dynamics of such a body be, if uniaxial stress is applied in addition, as in the setup of the 
Grinfeld instability? Of course, the assumption that the equilibrium stress remain constant is an 
oversimplification now, since the dielectric properties of the solid and its melt will usually differ, 
hence the electric field would become inhomogeneous as soon as an interfacial shape change occurs. 
Let us nevertheless assume the simplest possible situation, an anisotropic but constant equilibrium 
state 

o-^P^ = -poSij + XoSixSjx ■ (35) 

Using the stress-strain relationship (p7|), this can be inserted into our expression for the clastic 
energy density of the solid, which then becomes 

fsoiiiu.,}) = i^u.,u., + 2 u,, + Apuu + --^^^ xoun - Xo^^^ + Xo^j^^^^ , (36) 

where we have set poi = pq for simplicity. 

It is then straightforward to derive the sharp-interface limit for this modified phase-field model. 
The result reads (on setting ctoo — 0) 



- — < 7K -f Apgz + {att - cr„„ - xo) 



2 



^ Xo(o-it - o-«n) ^Xo 



+ ^ X^n^ } , (37) 



where is the component of the interface normal in x direction. Rotational invariance is broken. 

We are not aware of any previous mention of this equation in the literature, nor do we think this 
interesting case has been treated. In fact, what we have demonstrated here, is how the phase-field 
model can be used to derive hitherto unknown sharp-interface equations in a transparent way. 

It is clear from ( |3^ ) that an isotropic stress tensor, i.e., att = o'nn does not necessarily entail a 
stable planar interface, whereas setting att — o'nn = XO) i-e-, providing an anisotropic stress tensor, 
we will have a linearly stable planar front solution with interface position z — 0. This is easily seen 
from the fact that the terms containing and do not contribute in a linear stability analysis, 
because is directly proportional to the perturbation and hence its square and fourth power have 
to be dropped. Note also that the symmetry of the dynamics with respect to a replacement of 
Ctt — O'nn by its negative value does not hold anymore in this situation. 

While this equation opens a new line of research, we will refrain here from pursuing this topic 
any further. 



III. NUMERICAL RESULTS 



A. Validation of the model 



In order to verify that our phase-field description leads to a quantitatively correct description 
of the instability, at least before cusps set in, we have performed a number of numerical tests. 
Based on a simple finite-difference scheme, the numerical implementation is set up in a rectangular 
geometry. 

The bottom half of the rectangle is filled with solid, the top with liquid. This is realized by 
setting the phase field equal to a tanh-like function taking the value one in the bottom region 
and zero in the top region of the geometry, (j) is kept at these values one and zero exactly at the 
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bottom and top lines of the numerical grid, respectively. Periodic boundary conditions are applied 
at the lateral boundaries. The initial interface is set by an appropriate modulation of the region 
where (j) crosses the value i and was in most cases taken to be sinusoidal or flat with a random 
perturbation. 

The boundary and initial conditions for the fields and Uz are chosen differently for the KM 
and MG models as will be described now. 

Within the KM model, where we assume strains to vanish at equilibrium (hence Ap — 0), we took 
the X derivatives of both displacement fields periodic in the x direction in our initial simulations, 
in order to obtain periodic strains. Later, we switched to simpler helical boundary conditions 
for Ux, i.e., we took Ux{L,z) = Ux{0,z) + Luxxfi, where L is the length of the rectangle along 
the X direction, and periodic boundary conditions for Uz- This change in boundary conditions 
did not affect results in any essential way. All the simulations of the KM model discussed here 
were carried out with these boundary conditions (whereas those in were done with periodic x 
derivatives). At the bottom of the system, the values of the fields are fixed to values corresponding 
to a homogeneously strained solid; at the top, Ux is fixed and the derivative dzUz chosen such that 
the condition Vu = is satisfied. Ux is initialized as a linear function Ux = xUxx,o and the inital 
Uz is determined via integration of Eq. (A38). 

For simulations of the MG model (or rather its variant considered here) , the fields were all taken 
periodic in the x direction, whereas the boundary conditions in the z direction were as in the 
KM model. We did not yet attempt to use spectral methods for the solution which would require 
periodicity in the z direction as well (achievable by simply reflecting the system at its bottom, 
and including the image into the nu meri cal box p7|). Initialization was done by setting Ux = 
everywhere and computing Uz from ( [A3^ ) again. 

The elastic equations were solved by successive overrelaxation, the time integration was per- 
formed by a formally second-order accurate midpoint scheme. Since we did not update the elastic 
fields at the half time step, the formal accuracy was not attained. The most time-consuming 
part of the simulation was the relaxation scheme and a way to overcome its restrictions has been 
given in [ pT[ as is discussed in appendix B. Since it requires an approximation to the solution 
of the elastic problem even at the analytic level {n/K has to be small), we did not implement it 
in our two-dimensional simulations. We intend to compare the quality of this approximation to 
the solution of the full problem, before employing it in a 3D simulation, where its use is essential 
for reasons of computational efficiency. Most of our computations were done with the material 
parameters of Helium to facilitate comparison with experiments Therefore, whenever we do not 
indicate different choices, our parameters were chosen as described in [p6| . Times and lengths given 
without units are in seconds and centimetres, respectively. Since our nondimensional time unit is 
about one second and the nondimensional length scale about 0.1 cm, this simply corresponds to 
using 10^1 instead of ii as the basic length scale. 

One of our numerical tests consisted in reproducing the instability threshold to within 2% accu- 
racy, another one in verifying the subcritical nature of the bifurcation, first demonstrated analyti- 
cally by Nozieres ||2^. A short discussion of the last feature has been given in jl^], so we will not 
elaborate on it here ||29|]. We consider a few more tests, however. 

Figure 1 gives the dispersion relation determined for three values of the external stress and 
compares it with the analytic result from linear stability theory. The KM model was used here as 
it gave more accurate results at finite e. 

To obtain the dispersion relation, we simply followed the dynamics of a system initialized with a 
small- amplitude cosine profile for a number of different wavelengths, and computed the amplitudes 
of the evolving structure for a series of times. Then the amplitudes were fitted to an exponential 
function which provided the growth rate of the interface. Taking aoo equal to the applied stress 
ctqi we fixed the average position of the interface. Moreover, the amplitudes were computed in two 
different ways, both of which are not influenced by the average interface position. The first method 
was simply to take the square root of the spatial variance of ^ ^ second measure for the 

amplitude we took the modulus of the Fourier component corresponding to the wavelength chosen. 
On the figure, these two methods give essentially indistinguishable results within the size of the 
symbols. System sizes used were the wavelength A/ of the fastest-growing mode and a number 
of rational multiples and fractions thereof (ranging from jXf to 3A/). Since we kept the number 
of numerical grid points the same for all the systems at A/, the mesh size had to be varied. The 
interface thickness e was in general kept above | of the mesh size, which gives a resolution of 5 
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points for the region where the phase field varies between 10% and 90% of its maximum value. For 
smaller values of e, locking effects to be discussed shortly became conspicuous | |3l[ |. 

The agreement between analytic results and numerically determined points is satisfactory both 
above (ctq = 2.8 x 10^ dynes/cm^) and below (ctq = 2.4 x 10^ dynes/cm^) the instability threshold. 
Two points are worth mentioning. First at g sa 30 cm^^, there are two symbols each for the 
growth rates corresponding to the two larger stresses. These were given to roughly indicate the 
possible error in the numerical result when the growth rate has a large negative value. Points below 
g « 20 cm^^ did not show a comparable error. The two different values were obtained by fitting 
with the initial and the final half of the data points, respectively. We ascribe the difference to the 
fact that the amplitude of the interface becomes smaller than its width e, a situation in which the 
phase-field description is no longer reliable. For example, the final planar interface is not located 
exactly at z = 0, about which the initial cosine was centered, albeit the deviation is smaller than 
the interface width. Second, the overall agreement is surprisingly good in view of the fact that a 
phase-field model is not particularly well-suited to the determination of a dispersion relation at 
all. For in order to approach the limit of an infinitesimal perturbation of a planar interface one 
should choose very small amplitudes, but they must not be smaller than the interface width e. 
Reduction of e is possible in principle but soon leads to prohibitive computation times. With a 
sharp-interface model that we investigated in parallel | |30| |, it was no problem to take amplitudes 
of 10^"* and to obtain nice exponential growth or decay during long time intervals, whereas here 
we were restricted to starting amplitudes on the order of 0.05 or larger. 
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FIG. 1. Dispersion relation. Symbols indicate the results of numerical simulations, lines depict the 
analytic theory. Material parameters are those of Helium. Solid line and squares; ctq = 2. 8x 10* dynes/cm^, 
e = 0.009; dashed line and inverted triangles: ao = 2.6 x 10* dynes/cm'^, e = 0.01; dash-dotted line and 
triangles; ao = 2.4 x 10* dynes/cm^ e = 0.012. Mesh sizes; h = 0.0054, h = 0.0063, h = 0.0074, 
respectively. 

Our next test consists in investigating the dynamics of a planar interface with both the KM and 
MG models. From (|l|, |^) we obtain the equation of motion 

+ '''' 

which is, given the initial condition ^(0) — 0, solved by 

This analytic result is compared with simulations of the two models in Fig. 2. 
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What is cleared up by the figure is that even with a well-resolved interface width (we have 
e = 4/i) the MG model is slightly off the analytic final position, whereas the KM one converges 
well towards it. With larger values of the numerical mesh size, convergence of the former model 
gets even worse. For h = 0.007, e = 0.011 the KM model still agrees reasonably well with the 
analytic curve while the MG one is off by about 10% for t = 4. Both models show deviations 
from exponential behavior with this set of parameters due to metastability effects of the discrete 
set of interface points. This problem, which is particularly critical for interface pieces parallel to 
one of the coordinate axes, has been discussed in detail in At small interface velocities, the 
sum of the energies of the discrete points of the phase field in the double well potential may vary 
at successive time steps (whereas the energy of a continuous field is degenerate under arbitrary 
translations). Therefore, the interface is slowed down, if the energy of its discretization increases 
due to motion and accelerated if it decreases. For sufficiently small driving force, the interface 
may stop moving at all, i.e., lock into some favorable position. Apparently, the MG model is more 
susceptible to these effects than the KM one. 

The ultimate reason for the different behavior of the two models is that they are only asymp- 
totically equivalent, i.e., they describe the same system only in the limit e 0^. For any finite 
e, the equations obeyed by the phase field and the displacements are not the same in the two 
models. One can observe this directly by comparing the different terms contributing to, e.g., dt4>- 
In the MG model, the term Apuu of eq. (|2^) is frequently the largest interface term affecting dt(j), 
whereas this term is equal to zero in the KM model. Moreover, the sum of all terms multiplying 
h'{(j)) is not the same in both equations. 
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FIG. 2. Recession dynamics of a planar interface. Solid line: KM model; dashed line: MG model; 
dash-dotted line: analytic result (M). The dash-dotted line is hard to see; it almost coincides with the 



solid one. To discern it, one should look at the left part of the figure. 
h = 0.002. 



2 X 10* dynes/cm^, e = 0.008, 



The difference can also be seen in comparing a numerical simulation of the the sharp-interface 
model (|l|,|) itself, using an integral equation approach, with the phase-field models. We will report 
on details of this alternative approach elsewhere . Figure 3 shows the interface evolving in the 
phase-field calculation for two different values of e and compares them with the sharp-interface 
result starting from the same initial condition, after the same time interval. 

Again the KM model fares slightly better in the comparison for the same value of e. In the 
groove, however, both models deviate from the sharp- inteface result but approach it more closely 
for the smaller interface thickness. The sharp-interface model produces a more strongly pointed 
groove, as expected. It should be emphasized that this simulation is not far from the limit of the 
temporal validity range of the sharp-interface model. This limit is signalled either by a crash of 
the program due to the singular behavior of the bottom of the groove or by the appearance of a 
spurious steady state, which can be achieved by overstabilization of the interface. 
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Why the KM model agrees better with the sharp- interface model for a given value of the interface 
width is a difficult question, to which we cannot offer any deep answer. Also we cannot exclude 
that for a different choice of the functions h{(j)) and g{(j)), the MG model would be superior. It 
should be kept in mind that the functions employed in are not the same as those used here 
(see App. B). 

Normally, we would thus prefer to use the KM model for calculations to be presented. However, 
since we wish to make sure that effects of translational symmetry breaking are not due to our 
using a model in which periodicity is imperfectly implemented, we will use the MG variant in the 
following simulations. The differences between the two models are small, after all. Also the MG 
model has the advantage to be more easily treated using pseudospectral methods based on Fourier 
series due to its periodic boundary conditions, with a gain in accuracy that might allow to offset 
its apparent disadvantage. 




-0.08 - 
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FIG. 3. Comparison of phase-field models for different interface widths with the sharp-interface model. 
Solid lines: KM model; dashed lines: MG model; dash-dotted line: sharp-interface model. The phase-field 
curves with the shallower minima correspond to e = 0.01 (mesh size h = 0.0688), those with the deeper 
minima to e = 0.005 {h = 0.0344). ao = 2.5 x 10* dynes/cm^, t = 0.25. Initial interface amplitude: 0.05. 

The conclusion from Fig. 3 is that the phase-field models give decent agreement with a sharp- 
interface calculation in regions where the curvature is not too large. Whereas the sharp-interface 
computation cannot be meaningfully continued by very much beyond the time shown in the figure 
(t = 0.25), the phase-field models both have no problem in continuing the simulation to times well 
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beyond t = 1. 

As anticipated above, we take the point of view that a real sohd cannot develop exact cusps, 
because plastic effects such as the generation of dislocations will intervene. These will relieve 
stresses and thus prevent infinite densities of the elastic energy. The phase-field model does the 
same thing and we shall see below that it does so by introducing a cutoff to the curvature. More 
quantitative modeling would require to explicitly take into account models of nonlinear elasticity 
or plasticity, which is beyond the scope of this article. Nevertheless, as we can see from Fig. 3, 
the behavior far from the sharp tip of the groove is described reasonably well by the phase-field 
model for both values of e and is almost independent of the interface width. Therefore, we believe 
that the phase-field approach correctly reproduces the qualitative behavior of a situation in which 
plastic effects occur only in the minima of the grooves. 

Results obtained under this hypothesis will be discussed in the next section. 

B. Dynamics of extended systems 

When simulating periodic structures, one realizes that for small supercritical stresses, where the 
system takes a long time to develop deep grooves, one often observes symmetry breaking and one of 
the grooves getting ahead of the others. This symmetry breaking must be triggered by numerical 
noise from roundoff or tnmcation errors. For high stresses, where the system develops grooves 
reaching the system bottom within a relatively short lapse of time, this does not happen. Figure 
4 gives an example of a structure grown at about three times the critical stress. The interface is 
plotted at constant time increments {At = 0.05). A shift in the chemical potential of the solid has 
been made to keep the position of a planar front fixed. 
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8 X 10* dynes/cm^, e = 0.02. MG model with prestress o-qo = 8 x 10"* 
The fastest-growing mode is at A/ = 0.066, the wavelength of the pattern is A = |. Note that 



FIG. 4. Dynamics for ctq 
dynes/cm^ 
the two a:xes are the same scale 



We see that the structure remains periodic in the time interval considered and that three equally 
deep grooves evolve. Note the peculiar shape of the cells. Prom flat tips there emerge slightly curved 
slopes on the side of the cells. Then there is a sharp bend downward into the deep groove. The 
appearance of this bend renders it plausible that the time of formation of a cusp in the sharp- 
interface description has already passed and from then on the dynamics should be governed by the 
curvature bound. In the final stage of this dynamics all grooves move at constant velocity. Figure 
5 gives the curvatures of the interfaces displayed in Fig. 4 and demonstrates that the radius of 
curvature at the bottom of the grooves remains constant and is close to e, which was equal to 0.02 
in this simulation. 

The curvature was calculated from the contour line defining the interface position. Since the 
representation of this line {cj) = 0.5) was constructed by determining its intersection points with 
the squares of the numerical grid, the discretization points were unevenly spaced (two intersections 
with grid lines parallel to the x and y axes can be arbitrarily close to each other, the next may 
be as far away as V2h). Therefore, our curvature results arc pretty noisy, even after application 
of a smoothing procedure. A superior method for their determination would be to use the full 
representation or the phase field instead of just the contour line information. Nevertheless, they 
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clearly indicate the approximate constancy of the curvature in the groove tips. 




FIG. 5. Evolution of interface curvature, corresponding to the interfaces of Fig. 4. 



Since we have the stresses at our disposal, too, we can calculate the final velocity of the grooves. 
Figure 6 gives a contour plot of the stresses a^x , corresponding to the final time of the simulation 
from Fig. A, t = 0.45. The interface is drawn as as solid line, the contour lines are broken lines 
in different sty les. What we have plotted here, is not a generalized stress tensor component, as 
defined by (A4), but simply the stress in the solid. Therefore, the contour lines for stresses far in the 
liquid are meaningless (in the dynamic equations, they are multiplied by h{4>) ~ 0), although they 
become important when entering the interface region. From the figure, we estimate a maximum 
value of (Ttt « 2 X 10^ in the bottom of the groove (and a similar value is obtained from the 
corresponding figure for a^z)- Inserting this, the value of the curvature and the position ( of the 
groove bottom into (^), we obtain for the velocity w„ = —1.6. Assuming that the interface grew at 
this velocity from the outset, we obtain for its final position the value ( — —0.7. The data show 
that it is actually at ^ = —0.72, which is easily explicable by the inaccuracy of our estimate of 
the maximum stress. From the contour plot we do not obtain more than a rough figure as stresses 
vary rapidly in the interface region. 

In a straight and narrow crack, the stress scales with the square root of the distance from its 
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tip |2J|. Therefore, a reductfon of the tip radius by a factor of two will increase both the stress 
terra and the curvature term of (|^) by a factor of two as well. As long as the gravity term in 
is negligible (which, incidentally, it is not in the simulation of Fig. 4, its contribution is about as 
large as that of the curvature for the last curve), this means that the velocity of the groove will 
roughly double when e is halved. This trend has been confirmed in the simulations, although the 
observed ratio is slightly smaller than the predicted one, but then our grooves do not yet really 
have an extremely small width compared with their length. 




50 100 150 200 250 300 

x/h 

FIG. 6. Stresses ai t — 0.45 (lowest curve in Fig. 4). Six contour values are displayed, the sequence 
of line patterns is dashed, dash-dotted, dotted with increasing value of a^x and an increment of 4 x 10* 
dynes/cm^ between successive curves. To distinguish lines with the same pattern from each other, the 
values have been explicitly marked at those curves where there was enough room, e.g., for a^x = 4 x 10* 
dynes/cm^. The lower dash-dotted curve corresponds to the value of the applied external stress. 

The next three figures show a simulation at a stress roughly 20% above the critical value. Our 
numerical box contains six wavelengths of the pattern initially. One of the grooves has however 
been made by 2% deeper than its neighbours. Contrary to the situation in Fig. 4, no prestress was 
applied, so a planar interface would move downward to a new equilibrium position. This kind of 
motion is superimposed on the shape-changing dynamics and serves nicely in separating the curves 
on the plot. 
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FIG. 7. Dynamical evolution of a perturbed interface at ao = 1.16 CTc ((Jq = 3 x 10* dynes/cm"^, e = 0.025, 
E = 3.2 X 10* dynes/cm^). MG model, no prestress. Time interval between curves 0.25; the final time is 
2.5. After an initial phase of acceleration, the interface slows down and approaches a cycloid-like curve. 

The temporal dynamics can be divided into several stages. At first, the sinusoidal pattern 
changes its shape in the way already discussed by Nozieres |2^: the tips become flat, the grooves 
pointed. After some time, the interface becomes similar to a cycloid but with different depths of 
the grooves. Also, the dynamics almost comes to a halt. Below, we shall discuss the similarity with 
a cycloid in more detail (see Fig. 10). It holds up to t = 2.5 approximately, which is the time of 
the lowest curve in Fig. 7. At this point the apparent periodicity of the pattern has doubled. (Of 
course, strictly speaking this periodicity has been broken from the outset by our making one groove 
a little deeper. But this was only to avoid its being broken by numerical noise in an uncontrolled 
manner, i.e., to introduce a well-defined perturbation.) 

The groove that was ahead initially, wins the competition for the elastic field; the losing grooves 
fall back and even close again. This is shown in Fig. 8, displaying the temporal continuation of 
Fig. 7. In the initial structure of Fig. 8 (the solid line that is shallowest in the big grooves), the 
smaller grooves are deeper than in the final one (the dashed line which is deepest in the big grooves 
but shallowest in the small ones). At the end of the period of time depicted in Fig. 8, there are 
three clear survivors and three losers of the competition. 
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FIG. 8. Continuation of the evolution from Fig. 7. Initial time t = 2.75, time step At — 0.25, final time 
t = 5.0. The biggest groove is still accelerating, while the other even-numbered grooves have a roughly 
constant velocity and seem to decelerate towards t = 5.0. Odd-numbered grooves retract. They are deepest 
at t = 2.75 and have almost closed at t = 5.0. 

Finally, as shown in Fig. 9, only one groove survives. Its velocity is almost constant over a range 
of times. Eventually it slows down and grows sideways towards the end, which may have to do 
with the fact that it gets too close to the bottom of the numerical box (which is at ^ = — 1). Also 
gravity has a decisive decelerating effect here. 

What we observe, then, is a coarsening process that seems to proceed via imperfect period 
doubling transitions. Because our system has only six grooves, we cannot explicitly sec more than 
the first period doubling here. These transitions are local in the following sense. Not all grooves 
surviving the first period doubling get ahead of the others simultaneously. Rather what happens 
is that first the winning groove gets ahead of its nearest neighbours, screening each of them off 
the stress field on one side a little. This causes these neighbours to grow more slowly, making 
them screen off their next neighbours on the other side less. So these get ahead of their neighbour 
grooves, and so on. The perturbation made by one groove moves through the array in an alternating 
fashion. In an infinite system, one could imagine a series of "near" period doublings propagating 
through the array. These morphology changes are not exact period doubling transitions, because 
there is no restabilization of a structure with doubled periodicity. The system remains dynamic 
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(but see the discussion on gravity below) which means that the foremost groove does not get slower 
than its competitors, which would be necessary for length adjustment. 
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FIG. 9. Continuation of the evolution from Fig. 8. Initial time t — 5.25, time step Ai = 0.25, final time 
t = 8.0. The winning groove has a constant velocity most of the time. 



The first of these period doublings may be discussed analytically in some detail. Consider the 
shape of the interface close to the last time of Fig. 7. It can be modeled approximately by a curve 
that we would like to call a "double cycloid" . A parameter representation of this curve is given by 



X = ^ — Asinki, — B sin 2k^ 
z= — Acqs — B cos 2k$, 



(40) 
(41) 



Figure 10 compares a double cycloid with the interface aXt — 2. The wavenumber 2k {— 9.425) is 
given by the basic periodicity of the initial interface (before it is perturbed), the amplitudes A and 
B have been fitted "by eye" and the double cycloid has been shifted using translational invariance 
in the x direction. (Its position in the z direction can also be adjusted, which corresponds to a 
particular choice of the initial chemical potential of the solid.) 
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10. Comparison of the interface ai t — 2 from Fig. 7 (dashed line) with a double cycloid (solid 



Since we made only one of the grooves deeper than the others, the agreement of the groove 
minima is not quite perfect, as we can adjust only the depths of this groove and its nearest 
neighbours by an appropriate choice of the two constants A and B. Had we taken an initial 
perturbation of periodicity length 27r/fc instead of a local one in the simulation, a much better 
agreement would have been obtained. The purpose of this comparison, however, is not to claim 
that the interface shape goes precisely to a double cycloid but only to show that it may be well 
approximated by such a curve, which can be considered a cycloid (with amplitude B) modified by 
a small perturbation of twice its wavelength. In our fit shown in Fig. 10, we have A « B/IO. 

Our key observation is then that we can solve the sharp-interface elastic problem for a double 
cycloid exactly in an extension of the work of Gao et al. |24|, using a conformal mapping technique. 
This solution is given in appendix C. In what follows, we will neglect the gravity term, a procedure 
that we justify later. The evaluation of the nondimensional velocity via Eq. (34) for the double 



cycloid yields, in the bottoms of the grooves [see appendix, Eq. (C32)] 
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(42) 



where a = 'ik/qf is the ratio of the actual wavenumber of the basic cycloid and the wavenumber of 
the fastest-growing mode. The formula with odd m holds for the minima with depth —2Bk + Ak, 
that with even m for those with depth ~2Bk — Ak. A condition for the solution to hold is that 
there are no self-crossings of the curve, therefore we must require Ak + 2Bk < 1. Let us now assume 
that A <^ B, i.e. that the pattern actually is a slightly perturbed cycloid (where the perturbation 
has twice the basic wavelength w/k). Then the denominator in Eq. (^) in front of the braces 
goes to zero for even m as 2Bk approaches the value 1. This is the finite-time singularity, already 
identified by Gao al. |2j]. The velocity goes to —oo, if the braces remain positive, which they do 
for small enough a, i.e., when the wavelength is large enough. For small A, we can expand (^^. 
This gives 
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a formula that shows that the marginal value of a is 1. Thus, for wavelengths larger than that 
corresponding to the fastest-growing mode (a « 1), the velocity will diverge in the deepest min- 
ima, leading to cusps in the sharp-interface limit. We could leave the gravity term out of this 
consideration, because it never diverges for finite C. 

Now assuming we are at or slightly above the wavelength of the fastest-growing mode, we can see 



from (43) that for (1 — 2Bk)^ < Ak the velocity is positive in the secondary minima corresponding 
to odd m This means resolidification and closure of the corresponding grooves. 

Suppose for a moment that A — 0. Then the system with a sharp interface will evolve towards 
a cusped cycloid, i.e., 2Bk will increase towards 1. But this means that eventually a point will 
be reached where 1 — 2Bk is small enough that any perturbation will be larger than (1 — 2Bk)^ . 
In this case, our equations state that (for 1 — a ^ 1) the tip perturbed in this way in the right 
direction (i.e., the perturbation must reduce the depth of the groove) will recede again, its velocity 
will become positive. A groove tip that is perturbed in the other direction will approach the cusp 
singularity even faster and reduce the speed of its neighbours. Of course, not all perturbations are 
periodic; what happens when only a local perturbation is applied, can be seen from the simulation. 

What the analytic calculation shows, then, is that the first period-doubling bifurcation happens 
before the cusp singularity is reached, if the periodicity of the system is equal to the wavelength of 
the fastest-growing mode. Whereas the bifurcation to a set of alternatingly receding and advancing 
grooves may happen for any wavelength larger than this one, whether it happens before or after 
the predictable time of cusp formation will in general depend on the strength of the perturbations 
present in the system. In the simulation of Figs. 7-9, the periodicity of the unperturbed system is 
|, the wavelength of the fastest-growing mode is 0.5. 

Ordinarily, the time when the finite-time singularity appears in the sharp-interface system will 
be too short for the losing grooves to have appreciably retracted. In our phase-field model, there 
are no finite-time singularities, so the evolution can continue. It is then highly plausible that 
further period doublings occur, even though we have no analytic model for these. But on general 
grounds, we expect screening of neighbouring grooves to become more effective as all grooves get 
deeper. Hence the process should repeat, even at wavelengths far from, but above, A/. 

The difference between the cases of a wavelength close to that of the fastest-growing mode and 
one far above it is that in the former case, the first period doubling will happen before the time 
t* , at which cusps form in the sharp-interface limit, whereas in the latter case, it will happen 
afterwards. This case is, in fact, realized in Fig. 4, where the wavelength of the fastest-growing 
mode is about one tenth of the periodicity length. From Fig. 6, we can infer that the translational 
symmetry with respect to the basic wavelength A = | has already been broken by numerical noise 
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(the stress pattern does not show this periodicity in the upper half of the picture, this symmetry 
breaking will slowly propagate into the lower half where everything still appears periodic) . 

Another interesting conclusion from formula ( ^2[ ) is that for a > 1, i.e., for systems with small 
enough wavelength, stable steady states may be possible, because then surface tension may succeed 
in overwhelming the effects of stress. For a > 1 and A = 0, the formula predicts that a cycloid 
becomes stationary in its minima before the appearance of cusps. We hope to report on this aspect 
in the future. 

Finally, let us have a look at a system with a random initial condition. Figure 11 shows the 
evolution starting from an interface resulting from uniformly random perturbations of a planar 
front. We see that first some 10 waves develop, which is already a coarsened structure, as the 
wavenumber of the fastest-growing linear mode would correspond to about 24 waves fitting into 
the system. However, the initial amplitude is too small for this wavelength to become clearly 
visible. Some time later, there are much fewer features and eventually, only two grooves remain. 
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FIG. 11. Evolution of a random interface, ao — 3.5 x 10* dynes/cm^, e = 0.05. Periodic boundary 
conditions, i.e., MG model. <t < 11.6, At = 0.4. 



Whether one of the two will die off in the end is not clear, since this is a simulation with gravity. 
Hence the largest groove is bound to stop at some time, because the stress and curvature terms 



29 



remain constant once all other grooves are sufficiently small, but the gravity term continues to 
increase. If the second-largest groove still has a positive velocity when the first stops, it will not 
reverse its growth direction, but only grow to a point where its velocity becomes zero. 

An example, where the final state actually consists of two grooves, is shown in Fig. 12. Here the 
applied stress is smaller than in Fig. 11, so the pattern actually does come to a halt within the 
numerical box, after a long time (t « 60). Note that during most of the period where two grooves 
are dominant, one of them is ahead. Once it stops, the second approaches and in the end it has the 
same length as the first, to numerical accuracy. In the case of two periodically repeated grooves, 
this is to be expected for symmetry reasons. With three or more grooves, it is also conceivable 
that not all of them are the same length in the steady state. 
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FIG. 12. Evolution of a random interface with periodic boundary conditions, ao = 2.8 x 10^ dynes/cm^, 
e = 0.035, < t < 80.0, At = 0.8. The final interface lias only two grooves and no furtlier substructure. 
It is symmetric with respect to two symmetry axes at the two possible central positions between the two 
grooves. 

We think that in the absence of gravity, the situation in this strongly nonlinear region is very 
similar to the evolution of a Saffmann- Taylor finger in a Laplacian field. The Lame equations 
determining the displacement field are scale invariant just as the Laplace equation (and in fact, 
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Eq. ( |34| ) is scale invariant for £2 = 00). Once a strongly nonlinear state has been reached, none of 
the length scales discussed in section II can play a role anymore, since they only govern the local 
behavior of the growth pattern. The long-range elastic field will determine the factor au — CTrm of 
the destabilizing term in (34) and this factor will be the larger, the fewer competitors of a groove 
have grown to the same depth. This will lead to smaller grooves not growing anymore. This 
situation bears strong similarities to the growth of thermal cracks described in |Q. The main 
difference is that there a loser in the competition will simply stop growing. In our case, it will even 
shrink again, for the crystal can not only melt but also freeze again, and whether it will do so is 
simply determined by the chemical potential difference (Q). 

An analogous behavior is found in the side branching activity of a dendrite in the region about 



20 to 50 tip radii behind the tip JM 35|. There coarsening is observed, too, which also proceeds via 



imperfect period doubling. If this dynamics can be described in terms of a series of nonequilibrium 
phase transitions at all, these would have to be considered first-order transitions because of the 
discussed locality aspect. There is no diverging length scale in a single transition. 

We expect that the dynamics of large systems can be described by scaling laws similar to those 
given previously for the growth of needles in a Laplacian field |^^. The fact that "needles" can 
shrink again in the elastic problem should modify the long-time behavior of the needle density, 
which must pass through a maximum and then go to zero as a function of time, for any needle 
length. 

To some extent, this expectation is supported by the coarsening scenario described in [T^ ] in 
which extended systems without gravity are studied using random initial conditions. They mea- 
sured the Fourier transform of the height-height correlation function S{q,t) and observed dy- 
namical scaling. For early times, they observed a strong similarity between this behavior and 
early-stage soinodal decomposition in long-range systems. For later times, when the linear the- 
ory no longer describes the data, coarsening is evident: The location of the peak (/max of S{q,t) 
moves to smaller wavenumbers, as the peak height increases and sharpens. The peak height fol- 
lows S{qinax,t) ^ t""^-^, where a w 2, while the peak width sharpens with time as w ^ t~'' , where 
7 « 0.5. The former dependence is due to the total interface length increasing linearly with time for 
any unstable wavenumber. The latter dependence is due to competitive ordering between different 
wavenumbers, analogous to phase ordering. Within the accuracy of their study, they find that the 
structure factor shows scale invariance: S{q, i)/S'(qmax, t) = S*{q*), where the scaled wave number 
q* = {q- gmax)/w. Fitting to S* ^ {q*y and S* ~ (l/q*)'^, for small and large q* respectively, 
gives 6^1 — 2, and ip ^ 5 — 6. 

It is however difficult to assess to which time regimes these results correspond when compared 
with the present simulations, because the freedom to rescale parameters has been used extensively 
in [r^ . Since the vanishing of grooves does not seem to be a dominant mechanism of coarsening in 
their simulations, it is likely that the time windows considered in |17| and here have little overlap 
and that the stage of needle- like growth of the grooves is never reached in ||l^. Of course, the 
concept of short and long times is ambiguous in the absence of gravity due to the scale invariance 
of (^) for £2 = 00. However, one can compare the depths of grooves with their lateral distances to 
decide whether growth is best described as competition of wavenumbers (a Fourier space concept) 
or as competition of needles (a real-space concept). In this context, it is also important to realize 
that the validity of the simulations in is restricted to small values of 11 /K with the parameter 
r]o (see app. B) being 0(1), hence these simulations are quantitative only for external stresses 
(To = 0{ii). For these stresses, the state describable as a forest of needles is only obtained after 
a long time [of order (K/fj,)'^]. Therefore, the scaling exponents obtained in ||l^ might not be 
relevant to the scenario discussed here, which makes an analytic treatment along the lines of |Q 
even more desirable, because it could provide these exponents. 

In an infinite system what we have depicted here is probably just the continuation of the coars- 
ening scenario described in [ pT| . 

It should also be pointed out that for sufficiently wide systems, i.e., in particular for infinite 
ones, this dynamics may be an intermediate state only. Once a groove becomes sufficiently long, 
stresses along its side may become large enough to provoke a Grinfeld instability of the "side walls" 
of this crack-like structure, as has been shown by Brener and Marchenko [Q. Whether or not this 
happens, depends on how efficiently the stresses are relaxed along the grooves, on the speed of the 
grooves, on the perturbation amplitude, etc. This secondary instability might completely change 
the scaling behavior, possibly leading to tip splitting of the grooves and tree-like structures. So 
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far, we have not seen anything of this kind in our simulations. 

In a system of finite lateral extent (or periodicity) we think coarsening will in the absence of 
gravity generally lead to the disappearance of all grooves except one which will grow at constant 
velocity. If gravity is present, several grooves can survive and in a sufficiently deep system they 
will stop once they reach a depth where the gravity term compensates the stress one. 

IV. CONCLUSIONS 

In this article, we have constructed a class of phase-field models from a free-energy functional 
including the elastic energy density. A salient feature of the model is that the liquid is treated as a 
shear-free solid, which is to be contrasted with phase-field models taking into account hydrodynamic 
effects in solidification, where the solid is usually treated as a liquid of infinite viscosity . Our 
approach implies the artificial introduction of coherence conditions at the interface which is however 
counterbalanced by the fact that the only relevant elastic variable in the liquid is Vu. A whole 
class of models is obtained instead of a single one as a consequence of the freedom of choice for the 
state of reference used in measuring displacements. We compared the two most natural choices 
and found them to yield slightly different numerical results despite their asymptotic equivalence. 

Investigating a large number of laterally small and extended systems, we believe to be able 
now to describe the generic dynamic behavior. For systems smaller than the wavelength of the 
fastest-growing mode of linear stability theory but larger than that of the marginal mode (where 
surface tension stabilizes the planar interface), stable steady-state strucures are pos sible even in 
the sharp-interface limit. This is similar to the findings by Spencer and Meiron 123] for the case 
of transport via surface diffusion, even though we think the whole picture is more complex than 
what they described Here, we did not show detailed results on small systems but focused on 
extended systems. 

The case without qray ity is particularly simple, as the equation of motion can then be made 
parameter free [Eq. (^^ without the last term]. Initially, an interface may grow periodically but 
as soon as perturbations break the periodicity, coarsening will proceed via approximate period 
doubling transitions. Supposing randomized perturbations, the interface will, after a sufficient 
lapse of time, not look much different from one started with random initial conditions (compare 
Figs. 9 and II). If the system is of finite lateral extent (but infinitely deep), only a single groove 
will survive growing at constant velocity, determined by the final constant stress and constant 
surface tension terms near its tip. The final velocity will scale with the system width L and the 
radius of curvature e as u„ ~ L/e. All the other grooves will eventually retract, i.e., they will not 
even survive keeping a finite depth, which is different from the behavior of cracks ]|3^ ]. For wide 
systems, stresses near the groove tips may become large enough to trigger a secondary instability 
Pq] which would considerably modify the system behavior and allow the appearance of complex 
crack morphologies. It is however possible that this will arise only in the case of finite perturbations, 
as grooves may grow too fast in this situation for the instability to develop before it is "advected" 
(relative to the groove tip) into a region of very small stresses along the groove. If the system is 
laterally infinite, the system state will first follow dynamical scaling as studied in and should 
then cross over to the scaling dynamics described here with the number of grooves continually 
decreasing according to a power law, possibly with logarithmic corrections. Alternatively, the 
coarsening scenarios observed in ]p^ and in this article might be governed by the same scaling 
laws, with their difference being only apparent. The emphasis of was on the scaling laws 
governing coarsening, that of the present study is on the mechanism of coarsening. Obviously, 
this situation calls for large-scale simulations in order to determine the scaling exponents in cases 
where the mechanism presented here is definitely at work already. 

If the mentioned secondary instability becomes important, the identified state of competing 
grooves will be only of intermediate nature. Of course, all our considerations hold only as long as 
linear elasticity remains valid in the bulk, nonlinear elastic effects may alter the scenario. 

With gravity included (which was not considered in ]0), there are some modifications. First, 
it is now possible for a planar interface to be stable (apart from a vertical translation). Once the 
threshold of the instability has been exceeded, the behavior will be similar to the case without 
gravity. However, we predict that it is possible for several grooves to survive in a finite system 
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and that they will eventually stop growing, because the stress does not increase beyond a certain 
magnitude due to the lateral system width, whereas the gravity term increases as long as a groove 
gets deeper. That several grooves may survive has to do with the fact that now we have length 
scales in Eq. (|3^). More simple-mindedly we can immediately see that once the biggest groove 
stops, the second-largest will not retract, if it still has a downward velocity at that moment. Once 
the second stops, we can repeat the argument for the third, and so on. The final state will consist 
of a number of grooves, probably of different lengths and disordered. (For large-width systems, 
the aforementioned secondary instability may again complicate the picture.) In laterally infinite 
systems, it seems likely that a scaling state will prevail, possibly with a modified scaling exponent. 
Because now both stresses at the tips of the largest grooves and the gravity terms continue to 
grow, but they will both grow linearly with the length of the grooves. If initially the stress was 
large enough to overcompensate the gravity term, it will presumably stay like that. It is, however, 
not excluded that starting from specific initial conditions a system can be stabilized by gravity in 
the end. 
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APPENDIX A: DERIVATION OF THE SHARP-INTERFACE LIMIT 



In deriving the sharp- interface limit, we will restrict ourselves to the two-dimensional case as 
we did in discussing the sharp-interface equations (where we used only two stress components and 
only one curvature). The generalization to three dimensions is, however, straightforward. We may 
use ( p9| ) and ( ^2| ) as outer equations, to be used in the regions where the gradient of the phase 
field is small. For convenience, we set zq = 0- 

To obtain the inner equations, we transform to a local system of (orthogonal) curvilinear coor- 
dinates comoving with the interface, with one coordinate axis parallel to V0; the corresponding 
coordinate will be called r, while the second will be conveniently expressed by the arclength s 
along the interface |^0|. We introduce a stretched variable setting r — ep. It is then easy to see 
that a distinguished limit of (|29|), leading to a nontrivial inner equation that allows to satisfy the 
boundary conditions, is obtained by setting e — e. In saying this, we have assumed that the stresses 
and strains behave properly under rescaling, i.e., do not diverge. Designating by capital letters the 
values of the fields in the inner domain (where the gradient of the phase field is large) , we then 
have the inner equations 




25'(1>) 



4 



(Al) 

(A2) 
(A3) 



where 



= /l($)Sa;3 - (1 - H^))P6o.f, 



(A4) 



is the generalized stress tensor of the two-phase system. In order to obtain ( Al - A3 ) , we have used 
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V =^ndp + tds , (A5) 

1 2 



^2 



d'p + -dp + di . (A6) 



Derivatives such as dxUx can then be expressed in invariant form as ea;(ea;V)u, which leads to the 
following relations for the strain tensor components in the new coordinates: 

Upp = i dpUp , (A7) 
Uss = dsU, + nUp , (A8) 
Ups = Usp = i (^dsUp + i dpUs - kUs^ . (A9) 

The next step consists in solving the outer and inner equations via an asymptotic analysis that 
leads to a globally valid approximation for small interface thickness e which approaches the sharp- 
interface equations as e — > 0^. To this end we expand both outer and inner fields in powers of 
e: 



X, z, t) = (f)a{x, z, t) + e4>i{x, z, t) + ... , (AlO) 
'{^\x,z,t) + tuf^ 



Uij{x, z, t) — uf,\x, z, t) + tu]^Mx, z, t) + ... , (All 



and 



(j){x, z, t) = s, t) = $o(p, s, t) + e$i(p, s, t) + ... , (A12) 
u,j{x, z, t) = U,,{p, s, t) = Ulf{p, s, t) + eU^fip, s, i) + ... , (A13) 

where, due to the transformation properties of tensors, we can think of the subscripts i, j as 
running either over the values {x,z) or (r,s) and (p,s), respectively. Our basic field equations are, 
however, equations not for the strains but for the displacement fields. Thus, the expansion of the 
Uij{p, s, t) induces one for the displacement components: 

Ur - Up{p, s, t) = C/(0) (p, s, t) + eC/(i) (p, s, i) + ... , (A14) 
= U,{p,s,t) = C/f <) + £[/(!) + ... . (A15) 

Now the physical requirement that both Ur and Us remain finite in the limit e — s- 0^ allows us to 
conclude from (A7) and ( |A9| ) that neither C/p*^-* nor can depend on p, hence 

[/(") (A16) 
C/i") (A17) 

Furthermore, we have matching conditions for 1 <C p <C e^^ that can be obtained from the inner 
and outer expansions by equating equal powers of e (and taking into account that the variable r 
is itself e dependent): 

^o{p,s,t) (j)o(r,s,t)\r=±Q , p^±(xi, (A18) 
<^i{p,s,t) ^[(j)i{r,s,t) + pdr(j)o{r,s,t)]\r=±Q , p ^ ±oo , (A19) 

f^S(P'*'^)^^S(^'*'*)|r=±0, p^±oo, (A20) 

where we use the ~ symbol in the sense of asymptotic equality, i.e., f{x) ~ g{x), x — *■ Xq is 
equivalent to limj,^^,^ f{x)/g{x) = 1, and for two series in x — xq we require this relation for each 
corresponding pair of terms. 

The relations induced by ( A20 ) for the displacements are more complicated. We just give two 
examples. Because each derivative with respect to r comes with a factor 1/e when transformed 
into a derivative w.r.t. p, we have u'^p — dpU'^^ and hence 
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lim dpUj,^\p,s,t) ^ drul°\r,s,t)\r=±o . (A21) 
Our second example is even more instructive. We write 

lim C/(0)(p,s,i)= lim (dM^^Hs,t) + KUj,^\s,t)) 
= a.C/W(s,t) + «;(7(o)(s,i) 

= [dsuf^ [r, s,t) + Kuf^ (r, s, t)] |,=±o , (A22) 

which shows that the hnear combination u's^ + nu'p^ must be continuous across the interface. 
Finally, we need the expansions of all functions of cj) in powers of e, e.g.: 

h{cp) = /i(0o) + e/i'(0o)</'i + {h'{(t>o^2 + \h"{M(fi^ + ... , (A23) 

and we will use the obvious abbreviations /iq, /ig, etc., for functions of (pQ. Let us note a few useful 
relations in passing: 

(A24) 
(A25) 

)h"{P) , (A26) 





= 6(/)(l - 0) , 








= 20(1-0)(1- 



1/ 

18 



We have now collected all the prerequisites to perform the asymptotic analysis providing the 
sharp-interface limit. 

First, we note that the outer solution to lowest order [order e^^ of (p9|)] is simply given by 
g'{(j)o) = 0, which yields the solutions (j>o — 0, (f>o = 1, and (po = ^ [see ( A26 )]. The last of these 
is unstable and also not compatible with the boundary conditions in typical numerical setups. We 
assume 4>q = for r > 0, corresponding to the liquid phase, and 4>q = 1 for r < 0, corresponding 
to the solid phase. Equation (A24) tells us that h'{(j>o) — 0, and hence these solutions are valid at 
all orders of e. Using ft.((/)o) = in the liquid and /i(0o) = 1 in the solid, we immediately see that 
(|3^ ) turns into the mechanical equilibrium condition for the liquid and solid, respectively: 

dip = (hquid) , (A27) 
dja,j = (solid) . (A28) 

This is again true at all orders of e, and we can write the zeroth-order piece of the result in the 
form: 

= Poe - \ufl = po = const. , (A29) 
^ = -pos% + -Imf^ + \utl S,, . (A30) 

Later, we will look at two reference states in particular. One is the "natural" choice pos — pot, 
i.e., the unstrained state is hydrostatic and corresponds to the same pressure in the liquid and in 
the solid. If moreover, this pressure is chosen equal to the equilibrium pressure po, then we have 
uff! = in the liquid at equilibrium. The second choice corresponds to assuming a finite difference 
Ap ~ poi — pqs while keeping poi — pq. This means that zero strain corresponds to a prestressed 
solid, with a stress tensor cr.y = —poSij +ApSij, i.e., the deviation from equilibrium is the isotropic 

tensor Ap Sij . Both approaches can be exploited numerically^ 

We now consider the inner solution. The lowest order of ( Al ) gives 

dl<i>o - 2.g'($o) = . (A31) 

This equation can be solved by standard methods. Multiplying by 9p$0i we immediately obtain a 
first integral, written down here for further reference. 
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(A32) 



which is solved by <i>o — \ {\ — ta nh o), and this solution satisfies the matching conditions ( Al§| ). 
From the strain equations (A2), ( |A3| ) we obtain to lowest order 



5pS("; = . 

dptf^ = . 



(A33) 
(A34) 



which on integration from ootop^oo together with the matching conditions (A20) yields 

|.=o- = S W (-(X3) = (oo) = -P^ -PO , (A35) 
|.=o- = S W (-00) = (^) = . (A36) 

The limiting values of T,^pJ and T,^sp can be gathered from (A4). Obviously, these two equations 
constitute the condition of mechanical equilibrium at the interface, as dr? and air^ are the normal 
and shear components of the stress tensor of the outer solution there. 

However, the strain equations provide more information than just mechanical equilibrium on 
the outer scale. We write (A33) explicitly in terms of the strains and integrate indefinitely with 
respect to p, which yields 



(2^^ + A - A)C/(°) + (A - A)C/i°) + Apl + A(C/(°) + C/i°)) = f{s, t) 



(A37) 



where f{s,t) is a function of integration, to be determined from the matching conditions. This is 
straightforward and yields /(s, t) — poi — po- Moreover, we know that the only spatial dependence 
of Uss'^ is that on s [see (A22) and (A17)], which suggests to solve for C/pp-* in terms of [/is"*. The 
result is 



-1 



(2^ + A - X)ho + A 



Ap ho+po- Poi + [A + (A - \)ho]u!;°J 



(A38) 



The advantage of this equation is that it provides us with the full p dependence of C/pp' , allowing 
the explicit evaluation of integrals on p containing the strains. An analogous procedure for the 
second strain equation determines ui'P to be equal to zero. 

Now we proceed to the next-order equation for <i>. Written out explicitly, it reads 



Uf^) + Aiy + Ap.gz(s) 



(A39) 



where we have used = 0. With the help of ( |A24D and ( |A32D , we can arrange this as 



A- A 



(A40) 



with L = 9p — Ig'o being a self-adjoint linear operator. The solvability condition for this inhomoge- 
neous linear equation is that the right-hand side be orthogonal to the left-sided null eigenspace of 
h. Since L is hermitean, we know that the translational mode 9p<&o is an appropriate eigenvector: 



9p$oi = i5p$o = . 



(A41) 
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Multiplying ( A40) by from the left and integrating on p, we find 



dp\~kv + 7. +p[uif + uif] + + u^y 



+ Ap(C/(") + + AW- + Apgzis) \h',dp<fo . 



(A42) 



Now we can exploit (A38), telling us that the p dependence of the braces in (A42) is fully 
contained in their dependence on ho. All the integrals can be done analytically, using 



dpf{ho)h'odp<fo = - d^„fihi^o))h'i<i>o)= dhfih) 



Integrals that appear in (A42) are 

h = 
h = 
h = 



/QO 
dp /io<9p$o = 1 , 
-OO 

dpU^°^Kdp<fo, 
dpU^pfKdp^o 



(A43) 

(A44) 
(A45) 
(A46) 



The evaluation of the latter two integrals is as straightforward as that of the first, although a 
little more tedious. We just give the final result for the solvability condition, taking pQg = pq for 
simplicity: 



-kv = "fK + Apgz + 



2{fi + \){2p + \) 



2(/i + A)C/io)+Ap' 



2^ + A 
8p{p + A) 



(A47) 



At this point, we may specify our choice of reference state for the solid. First, let us assume that 
the unstrained state corresponds to a state of equal hydrodynamic pressure in the two phases, i.e. 
Pos = Po£i or Ap = 0. This is the KM choice. Then taking the limit p —f —oo of ( |A38| ) we get 



,(o)_4o), 



-Xu 



(0) 



2p + X' 



implying afs' — <Jrr' — 2p{u. 



(0) 



uiV) = 4/i(/i + A)ui?/(2/i + A), from which we obtain 



(0) _ 2/i + A (0) 

ss ^ A .. I \\^ tt "nn) 



Anin + A) 



(A48) 



(A49) 



where we have now switched to the conventional notation for the principal components of the 
stress tensor in the normal and tangential directions {<Trr = Cnn, Css — f^tt)- Finally, expressing 
the Lame constants by Young's modulus and the Poisson ratio, we arrive at 



1 n - 
I 



2E 



(4'-'^i°2)'-^oo + in + Apgz 



(A50) 



which is the desired sharp-interface limit. [In eqs. (|l|, ^), tioo = 0.] 

A remark is in order here. The phase-field equations imply the continuity of u'g} across the 
interface. As this quantity is obviously nonzero whenever the solid is strained, this means that we 
will not have u^ss = in the liquid. However, we will still have u~ss + u^rr = 0, i.e., the divergence 
of the displacement vector vanishes in the liquid. But this is all that matters, because it is only 
this quantity that enters the description of the liquid. 

The reason for ui} ^ in the liquid is that the phase-field description imposes coherence of the 
strain across the interface, which ultimately goes back to our viewing the liquid as a (shear-free, 
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but nonetheless) solid. For a true liquid in contact with a solid there is no such coherence condition 
as it is free to slip on the solid surface. Therefore, it could always keep its strain tensor isotropic 
(if such a notion made much sense at all for a liquid). The reason why we can nevertheless model 
the liquid as a solid is the additional degree of freedom that arises in the description of a liquid 
by having two fields Ux and Uz at our disposal even though only their combination Vu enters the 
free-energy expression. Therefore, we can compensate, so to speak, for imposing (nonphysical) 
coherence by allowing (equally nonphysical) anisotropic strain in the liquid. 

At this point wc may also note that if we were to model the elastic properties of two real solids 
by the current phase-field approach, we would necessarily impose coherence at the interface. The 
treatment of noncoherent solid-solid interfaces via phase fields would require some rethinking of 
the method. 

Concerning computational purposes, one disadvantage of the chosen reference state is that it 
is not very well-suited for the use of periodic boundary conditions, meaning periodically varying 
stresses and strains. The field equations are set up in terms of the displacements which acquire 
linearly increasing or decreasing components in directions where the strain has a nonzero average. 
This observation motivates the consideration of a different reference state in the solid, in which the 
average strain due to the external stress is subtracted. This is the MG choice (see App. B). If we 
impose a constant stress ctq in the x direction, the stress tensor in the solid is atj = — po^ij+fJo^ii^ij 
and requiring u^x = 0, we find that this is achieved by setting 



The corresponding homogeneous strain tensor is given by Uxx 
A). We then obtain taking the limit p — oo of ( A38) 



(A51) 

= 0, and Uzz = — Ap/(2^ + 



,(0) 



-Ap — Xu. 
2n + X 



(0) 



This can be used to express 



r(0) 



.(0) 



(0) 



2^1 + X 



2{fi + X)u 



(0) 



Ap 



(A52) 



(A53) 



wherefrom we obtain 2(/i 
leads back to ( A50). 



X)ufJ +Ap= {2^1 - 



A)(4°) 



o'nn)/2fi, which on insertion in (A47) 



Note that even here we cannot require u 



(0) 



^(0) 





at the interface and thus, according to ( |A38D , Urr = — Ap/(2(/x -|- A)), i.e., uiV would be constant 
along the interface. Then also cr^^^ — crin would have to be constant, which would lead to a dynamics 



in the liquid, which would imply u 

,(0) 



(0) 



entirely different from that of the Grinfeld instability [where (utV ~ ct„„ )^ increases in the grooves 
and diminishes on the peaks]. Hence, once again we are obliged to make use of the additional 
degree of freedom of the fields inside the liquid, even though now we can impose u'^x = Uzz — 
in the liquid as an initial condition (and as a far- field boundary condition), because this satisfies 
the periodicity requirement. 



APPENDIX B: MAPPING OF THE MULLER-GRANT MODEL TO THE PRESENT 

FORMULATION 

The main difference between the form of the MG model given in [ p^f and the one given here is 
a different choice of the functions g{<i)) and h{<^\ To make this conspicuous, we will rename their 
original functions to g{<i>) and /i(0) (since the second function was called g((\)) in [T^ , our renaming 
is also useful to avoid unnecessary confusion here). We shall leave the gravity and shift terms, 
/grav and /c out of the consideration, since they were not used by MG. 

Their double well potential is defined as 

/dw(<^) = -m , (Bi) 



38 



with §{(/)) — 0^(1 — (p)'^ , which is a sixth-order polynomial and actually has a third minimum at 
= — 1. The latter does not, however, play any role in the dynamics, provided no negative (j) 
values are given in the initial condition, a is a constant to be identified via the sharp-interface 
limit. 

Second, there is an elastic contribution to the free energy which they give as 



(B2) 



where K is the bulk modulus and fl the shear modulus which is (j) dependent: 

fj. = fiih{4>). 



The convenient choice 



(B3) 



(B4) 



guarantees that both bulk phases keep their equilibrium values at </) = (liquid) and (/) = 1 (solid) . 
This is due to the fact that h!{0) = h'{l) = 0, a property h{(j)) shares with h{ip) from the KM 
model (see Sec. II C). Obviously, the true shear modulus of the solid is /i = ^iA(l) 



For simplicity and since it does not change the behavior qualitatively, the bulk modulus is 
assumed to be the same in both phases. However, this restriction can be easily dropped by 
replacing K with 



K = Kn + Kih{(j)). 



(B5) 



As reference state they chose a prestressed state of the solid with Gxx = coj in which the strains 
Uxx and Uxz for a flat surface vanish. This entails that the state in which all strains vanish is a 
hydrostatic state with a different stress value. It is described by MG using a parameter 770 and 
[as may be verified easily from Eq. ( |B9| ) below], in this state we have (Txx = Ozz = ?7o^(l) = ^70/4. 
As a result, there is a relation between the new parameter 770 and the external stress cto in the 
uniaxially stressed reference state 



8(i^ + /ii/4) 

?7o = o-Q 

Ml 



(B6) 



The free energy density is then given as the sum of /dw(0) and /ci(0, {i*ii}) and two additional 
terms (77o/2iir)/i(0)^ and r]oh{(j))W ■ u. The first of these terms describes an energy shift, the second 
an additional coupling between the phase field and the elastic field (beyond that already implied 
by the dependence of the shear modulus and, possibly, the bulk modulus). A nice feature of the 
approach given in the present paper is that these terms are automatically generated by accounting 
for the fact that the equilibrium state does not have vanishing strain when the MG reference state 
is used: these are the terms containing po — pos in (|l^) and the corresponding terms Apuu and 
AW in (11). 

The free energy density is then given by: 



1 



2K 



h{d^f + 7]o h{(^)V ■ u + Ik{V ■ u) 



d 



(B7) 



The first term is the double well potential. The second and third terms are due to the partic- 
ular choice of reference frame, and 770 is related to the externally applied stress as described by 
Eq. (p36D. Applying the same line of reasoning as in Sec. II C, they obtain a system of coupled 
partial differential equations: 



-3'(0) 
a 



|/.(0)/7'(0) 



d 



(B8) 
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and 



dx, 



2^1 



d 



0. 



(B9) 



They show in that the phase field equations of this model also converge to the sharp interface 
equations. By expanding the solution of the mechanical equilibrium condition to first order in the 
shear modulus they were able to integrate out the elastic fields, so that they were left with an 
equation for only. That allowed them to use a pseudospectral method with which they could 
study wide periodic systems and three dimensional systems . 

Whether this expansion is entirely consistent is an open question, since they consider 770 an 
independent parameter that is 0(1), whereas for fixed ctq we actually have 770 = 0(l//^i), and /ii 
is the small quantity in their expansion. Probably this does not really matter in the absence of 
gravity where the precise value of 770 is immaterial as it only sets the time scale. The problem may 
be more awk ward in th e pr esence of gravity. 

Equations (B8) and (B9) are in a form that allows direct comparison with Eqs. ( p9|) and (^^, 
respectively. First note that in two dimensions ii' = A + /i = A + /ii/4 and that we must set \ — K 
to have the same elastic constants in the two sets of equations. In 2D, we have [u.^ — {5ij /d)W ■ u]^ = 



(summations over i and j are implied) and because of A — A = — /i, we obtain 



A- A , 1 



-Vu 



(BIO) 



which shows that on replacing h{(p) with jh((j)) the term of (BS) that is quadratic in the strains 
becomes equal to the corresponding term of (^^. Common prefactors will be discussed below. 
We then see immediately, that the choice Ap — rjo/i will make the linear term s equ al. This 
choice is the right one as is revealed by a quick comparison of Eq. (B6) with Eq. ( [A5l| ), derived 
in appendix A. Next we have to compare the constant terms which are rj^/K and AW . Here we 
notice that in |p^ , it was assumed that po£ = Po = 0. If we furthermore set ctoo = 0, then we 
infer from ^ that AW = Ap^ /2K = ril/i2K. Now the 0-dependent factor of -ql/ K in (|b|) is 
h{(j))h' {(j)) — ^d[h{(j)y]/d(j). We have checked by directly performing the sharp-interface limit of 
the original MG model, that this limit does not change when h{(j))'^ is replaced with |/i(0) (the 
solid and liquid phase limits are obviously unchanged). Doing this replacement first and then 
substituting jh{(j)) for h{(j)), we get identity of the constant terms, too. 

Finally, the prefactors should be discussed. In order to make the prefactor of the elastic ex- 
pressions the same in both equations, we must set F — l/(3fce), which shows that P = 3je. To 
determine the factor 1/a of the double well potential in (B8), one must actually perform the sharp- 
interface limit [because the potential is not the same as in (p9[)], which yields a — 3e/87. With 
these choices, Eqs. (B8) and (p9[) are asymptotically equivalent. 



The comparison of the equations describing mechanical equilibrium is even more straightforward. 
Inserting the expressions (|2^) and ( p8| ) with poi = and pos = — Ap = —rjo/A into (^) and using 

A ~ K, we get 



dx^ 
d 

dxi 



^% + KukkSij + 2^ (uij - -Ukk^ij 



+ [1 - h{ct>)]KukkS 



h{(j)) + Kukk 



2fj, 



d 

dxi 



hi<t>) u 



1 



UkkSv 



(Bll) 



which is obviously identical to (B9), once we replace h{(j))/'i by h{(j)). 



APPENDIX C: ANALYTIC SOLUTION OF THE ELASTIC PROBLEM FOR THE 

DOUBLE CYCLOID 

In this appendix (only!), we switch back from our notation for geometric relations in the plane 
which was based on a coordinate system spanned by the x and z axes (thus reminding ourselves 
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that in reality we have a three-dimensional system and we simply suppress deformations in the 
y direction) to the more conventional use of x and y for the planar coordinates. This way we 
"liberate" the symbol z for use as a complex variable: z = x + iy. 

We wish to solve the elastic problem in a half-infinite geometry, the top of which is bounded 
by the double cycloid given by Eq. (41) with z replaced by y. In the complex plane, this curve is 
described by 



z = x + iy = iAe ''^^ - iBe 



(CI) 



and because the parameter ^ is real, this equation also defines a conformal mapping frome the z 
plane to the C, plane, where C = C + mapping the curve z{x) to the ^ axis. 

To rephrase the elastic problem in the complex plane, we use the Goursat function formalism. 
Its basic statements can be easily inferred from the representation of two-dimensional elasticity 
in terms of a single scalar function, the Airy function x(a;,y). Setting Gxx = dyX, Oyy = d^Xi 
and Uxy = —dxdyX, the mechanical equilibrium equations djaij = are automatically satisfied. 
Hooke's law for isotropic bodies then implies that x obeys the biharmonic equation: 



V\ = 0. 

In terms of the complex variables z and z = x — iy, the Laplacian becomes 
dx = dz + dzj dy = i[dz — 9j)], hence the most general form of the solution to (| 

X{x,y) = x{z,z) = zfi{z) +gi{z) + zf2{z) + g2{z) , 



(C2) 

idgdz [because 
is given by 

(C3) 



where fj, gj {j = 1,2) are analytic functions of their arguments. Since we are looking for real 
solutions, we can restrict ourselves to two independent complex functions instead of four, i.e., we 
can write 



z (j){z) + / ip{z)dz + z (f>(z) + / il){z)dz 



(C4) 



where an overbar denotes complex conjugation. In this formula. 
From (C4), we obtain by direct differentiation 



i and if) are the Goursat functions. 



= V^x = 2 



'vv 



(/.'(z) + 0'(z) 



'vv 



+ 2iaxy^2[z<j)"iz) + ^P'{z)] 



(C5) 
(C6) 



The displacements can also be expressed by the Goursat functions [using Hooke's law and 
but since we do not need the corresponding relation, we omit it here. 

We have boundary condi tions fo r the stresses on the curve given by Eq. (CI), which one could 
attempt to use directly in (C5 C6) to obtain equations for the Goursat functions. However, since 
we are going to employ conformal mapping, it is useful to keep the order of derivatives small - 
the analytic evaluation of higher-order derivatives can be quite cumbersome. Thus it is desirable 
to reformulate the boundary conditions in partially integrated form. The force (fx-ify) on the 
boundary can be condensed into a single complex number 



fx ^" ^fy ^xj^j ^ ^^yj^j {^xx ^^xy)^x i^yy ^^xy)^^y j 



(C7) 



where (ux, Uy) is the normal vector to the boundary. Introducing the arclength s, we have (directing 
s such that s — s- cx) corresponds to x ^ cx)) 



and thus 



fx ify 



1 I 

0'(z) 

. d 

ds 



<j,'{z) 



ds 

. dz 

ds 
. dz 

ds 



. dx 



.dz 



ds ds 



i<^y 



(C8) 



z(j)"{z) 

[^(z) + z'^ + . 



4>'{z) 



. .dz 

ZKJxy ) 1 ~ 

ds 

dz 



ds 



(C9) 
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Integrating this local relation along the boundary, we obtain 



-z/ =-tJ (/, + tfy) ds = 0(z) + z0'(z) + i^iz) , (CIO) 

which allows to apply the boundary condition in expressions involving first-order derivatives only. 

We substract the stress at infinity from our (linear) elastic equations to be able to work with 
analytic functions that are bounded at infinity. Hence, we set 

'^v = '^ij^ + '^oSixSjx (Cll) 

and replace with ct^^' in Eqs. (C5,C^) above. Taking the equilibrium pressure in the liquid equal 
to zero (which we can do without loss of generality), the boundary conditions at the liquid-solid 
interface {(JijUj — 0) become 

translating into 

fx + ify^cFo^ -if =]^{z- z)aQ, (C13) 

where we have dropped an arbitrary constant of integration. 

Before embarking on the actual calculation, we ought to ponder one more point. We would like 
to solve a problem with periodic boundary conditions for strains and stresses in the x direction. 
But periodicity of the strains does not imply periodicity of the displacements nor does it imply 
periodicity of the Goursat functions. On the other hand, the use of periodic functions greatly 
facilitates the derivation. As noted by Spencer and Meiron |2^, we can express the Goursat 
functions by periodic functions 0o a-nd ■00 via 

(j){z) = (j)o{z) , 

^{z) ^ Mz) - zcjy'oiz) . (C14) 

The mathematical problem is then to find two periodic functions 00 and ip^, analytic in the 
domain occupied by the solid, satisfying 

0o(z) + (z - z)0o(z) + V'o(2) = -^(z - z)(To (C15) 

at the interface and remaining bounded for — s- — oo. The solution to this problem must be unique 
apart from possible additive constants to the functions (I)q{z) and "i/Wf )■ 



We transform to the C, plane, using the analytic continuation of (CI) 



This maps the interface to the real axis and the solid to the half plane rj <Q. To designate functions 
in the C, plane, we put a tilde on the letter they have in the z plane: 

Mz) - M^iO) ^ MO , 

Mz)^M<:)- (C17) 

The derivative of (po transforms as follows: 

Our task therefore is to construct two analytic functions satisfying 

MO + HO-m]=§r+l'oiO = -^HO-^m (C19) 
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on the real axis (77 = 0) and remaining bounded as — *■ —00. uj{^) is given by Eq. (C16), hence 

cj'(0 = 1 - Ake-'''^ - 2Bke-^'''^ (C20) 
and Eq. ( C19D becomes {A and B are real) 



MO 



(C21) 



Note that 0o(^) and ■0o(O ^'"S functions that should be analytically continued to the upper half 

plane, for if (/)q(C) [V'olC)] is analytic in the upper half plane, 0o(C) [■*/'o(C)] is analytic in the lower 
half plan e wh ich is what we need. The basic idea in constructing the solution is to divide the 
terms in (C21) into two groups, one of which corresponds to functions analytic in the upper half 
plane, the other to those analytic in the lower half plane. The equality between these two groups 
implies that each of them is equal to a constant, which gives us two equations for the two functions 
sought. It is clear that (^o(C) belongs to the terms of (C21) that are analytic in the lower half plane 

(on replacement of ^ by C,) whereas V'o(C) has to be analytic in the upper half plane. The difficult 
term is the middle one on the left hand side as it contains some expressions that are analytic and 
bounded in the upper half plane (e.g., e''^^) but also some for which this is the case in the lower 
half plane (e.g., e^*'^^). One way to proceed is to expand both o an d t/jq in a series in powers of 
e~'^^^ (a one-sided Fourier series, so to speak), to multiply Eq. (C21) by the denominator of the 
middle expression, and to separate terms with plus signs and minus signs in the exponents. This 
gives a two-termed recursion for 0o containing two constants that have to be determined from 
the analyticity properties. As it turns out, the series for 00 is finite, only the two first terms are 



nonzero. This suggests that a close look at Eq. ( C21 ) would have revealed this property, allowing 
to avoid the tedious expansion procedure. Indeed, there is a more elegant way leading to this 
result. Its discovery is left as an exercise to the astute reader. Here, we simply take 



0o(O = aie 



026 



(C22) 



as an ansatz. Inserting this into (C21), we get 



aie "-^^ + a2 



^ ' 1- Ake^^^ - 2Bke^'''i 



~M0 



(C23) 



Since we have just il'oiS,) on the right-hand side, which must be analytically continued to the upper 
half plane and remain bounded there, the left-hand side must not contain, after substitution of 
for ^, any "dangerous" terms diverging as 77 — *■ 00. Dangerous terms would obviously be the terms 
e~'^'"' and e"^**^^ as well as the zeros of the denominator. Now, we have the condition Ak + 2Bk < 1 
and we have |e*'^^| < 1 for 77 > 0. Therefore, the denominator is always different from zero in the 
upper half plane. All we have to do then is to choose ai and 02 such that the dangerous terms 
cancel. This is straightforward for 02, since there are only two terms containing e"^''''', after 
the prefactor of the second term has been multiplied with the numerator. There remains a term 
proportional to e~'^'"' , however, in the numerator, which must also be canceled by the choice of ai. 
The result is 



ai 



2 1 - Bfc ' 



a2 = —iB . 



From this, we immediately get (j)o as 



A 



-ikC 



Bk 



Be 



~2ifcC 



(C24) 



(C25) 
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Once (f>Q is known, ipo is obtained from (C21) or, using (C24), from (C23). We give the result for 
reference purposes, even though it will not be needed in the following 



1 



k{2B^ + A- 



l + Bk 
1-Bk 



2ABk 
l-Bk^ 



1 -Bk 



(C26) 



Note that only exponentials of negative multiples of ikC, appear; the numerators are thus analytic 
and bounded in the lower half-plane and the denominator remains nonzero because of >lfc+2_Bfc < 1, 
hence ■00 (C) satisfies all the required analyticity and boundedness conditions. In the case ^ = 
or B = 0, these results reduce to those of Gao et al. for the simple cycloid (after transforming 
back to the nonperiodic Goursat functions). 

We need only (j}Q to compute the tangential stresses on the boundary. Since we know the normal 
stresses to be equal to zero from the boundary condition, we can write att = tra = a^J + ay^J + ctq. 
We have 



_^'(C) c.'(C) 

which, when specialized to the boundary, gives us 



CTO 



1 + i+||Afce-'''''« + 2Bke~^'''i 
1 - Ake-'^'i - 2J5fce-2*fe« 



+ c.c. 



(C27) 



(C28) 



In the grooves of the pattern, we have fc^ = tott hence e = ( — 1)™, e ^'^^^ = 1, from which we 
obtain the tangential stress as 



1 I l+Bk 



^Afc(-l)™ + 2Bk 



1 - Afc(-l)'" - 2Bk 



(C29) 



To obtain the normal velocity in the grooves, we need the curvature there. The curvature is 
easily calculated from the parametric representation of the double cycloid 



( fc2 + 4B2 fc2 + 4ABk^ cos k£, + l- 2Ak cos k^ - ABk cos 2k£,f^^ 
and its value in the bottom of the grooves is (cos fci^ — ±1) 

fc2 [A(-l)™ +45] 



[1 - Ak{-1)"' - 2BkY 



(C30) 



(C31) 



Note that as a cusp is approached {Ak + 2Bk 1), crfj and the curvature diverge with the same 
denominator, for even m. 

Inserting ( C29| ) and (C31) into Eq. ( p4[ ) for the nondimensional normal velocity and neglecting 
the gravity term we obtain 



1 



2[1- Afc(-l)™-2B/c] 



l + 2Bk + ^^^Ak{-iy 
1 — Bk 



2£ifc2 [^(-1)" + 4B] 



(C32) 



where we have nondimensionalized the curvature via multiplication by £i. Setting a — 2£ik, we 
generate Eq. (E|). Note that if we consider the amplitude ^ to be a small perturbation of a cycloid 
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determined by B, the basic wave number is 2k, not k. Then a is just the ratio of the basic wave 
number and the wave number of the fastest-growing mode. 
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